//Copyright 2010 Neurosciences Research Foundation, Incorporated
//make sure all the nodes have the following files:
// connect.dat
// hosts
// Make sure the subfolder /voltages/ and /bold/ exist
// Change NP - the number of processors
// Change Br

#include "backtrace.h"

#include <limits>
#include <unistd.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h> 
#include <iostream>
#include <map>
#include <string>
#include <algorithm>
#include "constants.h"
#include "alex_mental_rotation_world.h"	// the interact function to communicate with the robot.
#include "neuronspecs.h"
#include "microcircuitry.h"
#include "groups.h"
#include "neurons.h"
#include "kbhit.h"

// GSL random generator - synchronizing between Mac and Linux systems
#include <gsl/gsl_rng.h>

//#include "world_sim.h"

using namespace std;

gsl_rng * rand_generator; // GSL
const gsl_rng_type * randgen_type = gsl_rng_default; // GSL

#define nomad

#ifdef SIMULATED_ROBOT
#define EXP2 
#else
#define EXP3
#endif // SIMULATED_ROBOT


// To restart from a checkpoint, comment out the start_dynamics #define.
#define new_branches
#define new_somas
#define new_synapses

//#define start_dynamics
//#define Save_Dynamics_Every_Second

//#define LOG_GROUP
#ifdef LOG_GROUP
  const int LOG_COUNT = 1;
  int GROUP_AREAS[LOG_COUNT] = {1};
  int GROUP_TYPES[LOG_COUNT] = {1};
  int LOG_START_MS[LOG_COUNT] = {0};
  int LOG_STOP_MS[LOG_COUNT] = {128000};
  /*
    Example:
    const int LOG_COUNT = 2;
    int GROUP_AREAS[LOG_COUNT] = {3,3};
    int GROUP_TYPES[LOG_COUNT] = {1,3};
    int LOG_START_MS[LOG_COUNT] = {100000,100000};
    int LOG_STOP_MS[LOG_COUNT] = {102000,102000};
  */
#endif

#define Allow_Plasticity

#define WRAPAROUND

//#define WARNINGS


// MPI transmission options: ONLY ONE MUST BE DEFINED AT ALL TIME!!!
// Theoretically variable transmission is fastest, but as the network 
// changes it may be worthwhile to test this empirically once in a while
//
// FIXED = allgather, as it has been in the past in this code, fixed xmit of int indices every ms.
// VARIABLE = allgatherv, as it used to be in TC code, different xmit size every ms
// BITPACKED = as it is in PaSpNet, 1 bit per neuron every ms
//#define FIXED_TRANSMISSION
#define VARIABLE_TRANSMISSION
//#define BITPACKED_TRANSMISSION

// transmission specific variables
#ifdef VARIABLE_TRANSMISSION
int neuron_count;
int displacements[NP];
int spike_count_in[NP];
int *spike_buffer_in;
int *spike_buffer_out;
#endif
#ifdef BITPACKED_TRANSMISSION
int nbytes_packed_out;
int nbytes_packed_in;
unsigned char *byte_buffer_out;
unsigned char *byte_buffer_in;
bool *index_buffer;
#endif


// Write the middle cell from each group instead of all cells.
//#define WRITEVOLTAGES1CELL
// Define WRITEVOLTAGES if you want to write voltage data at all.
//#define WRITEVOLTAGES


#define MINIs (0) 		// the frequency of minis, 1.0 means 1 mini per second per synapse
const int homeostasis_update_seconds = 1000000;
const int homeostatic_scaling_adjustment_start_sec = 1000000;
const int save_synapses_period_seconds = 1000000;
const int voltage_writing_frequency_seconds = 1;



// BCM
/*
 const float Aminus=.1;
 const float Aplus=.2;
 const float tau_minus = 34;
 const float tau_plus = 14;
 */

// original model:
//const float Aminus=.2;
//const float Aplus=.1;
//const float tau_minus = 19.5;
//const float tau_plus = 19.5;
//const float Aminus=0.0;
const float Aminus=.001;
const float Aplus=.005;    // .2
const float tau_minus = 19.5;
const float tau_plus = 19.5;

// ST-STDP constants
const float Aminus_sds = 0.12;
const float Aplus_sds = 0.1;
const float tau_minus_sds = 20.0;
const float tau_plus_sds = 20.0;


//int		 	sim_seconds = 30;  // Set this in world.h				
//const	int		NP	= 8;					// The total number of different subnetworks (each per processor)

#define Stimulate  	// size 2*250,000, DELAY 20; 2200 sec - trantient, 1100 sec -stimulation
					// size 400,000, DELAY 20; incoherent (random) stimulation
					//#define Lesion			// size 3*250,000, DELAY 10;

//#define fMRI				// 


#ifdef Lesion
#define FindGlobalAxons
#endif
#ifdef fMRI
#define FindGlobalAxons
#endif




#ifdef nomad
#include "mpi.h"
#include <sys/time.h>
#include <unistd.h>

// Define the global communicator to allow all procs in this simulator to talk.
// Not using MPI_WORLD_COMM any longer so we can communicate with other simulators.
MPI_Comm network_comm;
MPI_Comm exchange_comm;  // Allow communication with external simulator running with MPI.

#ifdef DISPLAY_WINDOW
MPI_Comm display_comm;
const int display_max_transmit = 1024;
int N_spikes_for_display;
int sp_buffer[NMAX];
#endif

//#define getrandom(max1) (((rand())%(max1))) // random integer between 0 and max-1
#define getrandom(max1) (gsl_rng_uniform(rand_generator) * max1) // GSL

#else

#include <iostream.h>
#define getrandom(max1) (((rand()*RAND_MAX+rand())%(max1))) // random integer between 0 and max-1

#endif 


//#define random11() ( (1-2*((float)rand())/RAND_MAX ) )	// random real in [-1 1]
//#define random01() ( (((float)rand())/RAND_MAX ) )		// random real in [0 1]
#define random11() ( (float)(1 - (2*gsl_rng_uniform(rand_generator))) ) // GSL
#define random01() ( (float)(gsl_rng_uniform(rand_generator)) ) // GSL

int N;

#define Rc 1.0f							// Size (radius) of the cortical sphere
const	float	Rt=1.0;						
// Scale factor for thalamic size.  Assuming 3mm barrel ctx and 600 micron barreloids on the long axis (The Thalamus, v2).
// So the axon size is reduced 
const float cortex_to_thalamus_size_ratio = 3.0/.6;  // Approximate based on The Thalamus, v2.

int firing_rate[SIZE]={};				// spike counts by neuron used for homeostasis.  Reset each homeostasis_update_seconds second.

int exchange_rank;
int display_rank;

int myid = 0;					// The id of this subnetwork (to be determined dynamically)
char processor_name[MPI_MAX_PROCESSOR_NAME];

const	int		max_transmit=NMAX/NP/2;		// The maximal number of spikes transmitted by a single subnetwork to other subnetworks

int	N_spikes_out;		// The number of fired neurons at a given time t (ms) from this processor
int	spikes_out[max_transmit];	// The indeces of fired neurons at a given time t (ms) from this processor

int	spikes_in[NP*max_transmit];// all global spiking information. Available to everyone.
int	spikes_indeces_in[NP];// tells where to look for spikes addressed to a processor

const	int 	training_dim=3;

float 	training_array[100000][training_dim];  // input array of training data : e.g. force_dv -- read from function read_training_array

//  Define the object controlling interaction with the outside world.
//world motor_world;

#ifdef EXP3

class g_motor_world
{
private:
	g_motor_world()
	{
		
	};
	
public:
	static world* Instance()
	{
		static world instance_;
		return &instance_;
	}
};

#endif

#ifdef EXP2
world motor_world;
#endif

//  Allocation of variables and parameters of the network.

const	int		M=1500;						// The number of synapses per compartment.  
const	int		D=DELAY;					// The maximal conduction delay in the network


const	float	myelinated_speed = 2*Rc/(D/2); 

int		sec, t, simtime;				// global variables telling the simulation time=1000*sec+t;
int		start_sec = 0;				// global var for first second used for restarting sim from checkpoint.

int		B_firings;							// the beginning of the que with spikes
int		N_firings;							// the end of the que with spikes
const int N_firings_max=D*NMAX/2;					// upper limit on the number of fired neurons per sec
int		firings_index[N_firings_max];			// indeces and timings of spikes
short	firings_time[N_firings_max];			// indeces and timings of spikes

#define sp_memory 5							// the number of last spikes stored per neuron
int		somatic_spikes[NMAX][sp_memory];		// the timing of the last 'sp_memory" spikes per each branch
bool		spiked[NMAX];					// spiked[n] == true if neuron n spiked this time step.
float		exc_input[NMAX];				// exc_input[n] is current added to neuron n from world.
int		dendritic_spikes[Br][sp_memory];	// the timing of the last 'sp_memory" spikes per each branch
#define LT_plast 1000						// The lendth of LTP or LTD window of interest (must be > 50+D ms)
float   LTD[LT_plast];						// LTD part of the STDP function
float   LTP[LT_plast];						// LTP part of the STDP function
float	LTD_sds[LT_plast];					// eligibility trace for ST-STDP
float	LTP_sds[LT_plast];					// eligibility trace for ST-STDP



// Compartments
// Branches
int		N_br	= 0;				// The actual number of branches (used as index for initialization)

float	Br_I_injected[Br];				// AMPA conductance in the compartment

float	Br_AMPA[Br];				// AMPA conductance in the compartment
float	Br_NMDA[Br];				// NMDA conductance in the compartment
float	Br_NMDAVI[Br];				// NMDA conductance in the compartment
float	Br_GABAa[Br];				// GABA_A conductance in the compartment
float	Br_GABAb[Br];				// GABA_B conductance in the compartment
float	Br_GABAc[Br];
float	Br_V[Br];					// The voltage potential in the compartment
float	Br_U[Br];					// The recovery variable in the compartment
float	Br_Isyn[Br];				// DEBUG; This is for writing to voltages to debug network.

int	Br_N[Br];					// The mother neuron of the compartment
unsigned char	Br_L[Br];			// The layer in which this compartment is located
									// Links: -1 means no link
int	Br_Down[Br];				// The link to the mother branch
int	Br_Up[Br];					// The link to the daughter branch
int	Br_Sister[Br];				// The link to the sister branch


// TODO: do not store locations, but store only the distinct ones and then re-use them
// Locations
//const	int	Max_N_locations=250000;	// The max size of the file that stores locations 'connect.dat'
//		int	N_locations=0;				// The actual size of the file 'connect.dat'
//		float	xxx[Max_N_locations][3];	// locations from the file (coordinates of the some on the sphere)


// Neurons	
int		Soma[NMAX];					// The link to the soma branch (c,-1,u,-1)
int		Inputs[NMAX/2];				// List of neurons to stimulate.
int		N_Inputs;				// Number of neurons to stimulate.

unsigned char	Type[NMAX];			// the type of the neuron
unsigned char	Area[NMAX];			// the area the neuron is in.
float		x[NMAX][3];					// location (coordinates) of the soma on the sphere
										//		float   	*y;						// location of the axonal target in DTI data ( y[N] will be allocated and freed dynamically via malloc)
										//		short   	xydist[NMAX];				// the fiber distance between x and y, according to DTI data
										//		float		Gx[NMAX][3];				// location (coordinates) of the global axonal ending (related but not equal to y)
										//		int		Gpre[NMAX];					// The index of the neuron with the global axonal ending
										//		int		N_G=0;					// The number of global axonal endings 

float		STSP_s[NMAX][2];			// Two types of short-term synaptic plasticity strength variables  
										// Each neuron may have up to 2 variables describing the state of 
										// two types of synapses, e.g., RS->FS or RS->LTS. Which variable 
										// corresponds to which synapse is given by STSP_x_var vector defined below.
float	Background_AMPA[NMAX];		// Background excitatory conductance
float	Background_GABA[NMAX];		// Background inhibitory conductance injected into all compartments of a neuron for awake state.


// Synapses
short	N_pre[Br];					// The number of presynapses. typically M or M-1

//int	pre[Br][M];					
//unsigned char	pre_delay[Br][M];	
//float	*s_pre[Br][M];				
//float	*sd_pre[Br][M];			

int **pre;							// Indixes of the presynaptic neurons
unsigned char **pre_delay;			// Delay from pre neuron to this branch
float ***s_pre;						// Pointer to the synaptic strength
float ***sd_pre;					// Pointer to the derivative of the synaptic strength
float ***sds_pre;					// Pointer to the ST-STDP strength

// N+1 - brainstem input, N+2 - sensory input
short	delays[NMAX+2][D];			// Distribution of conduction delays from pre neuron to all post neurons
int	*post[NMAX+2][D];				// indexes of the postsynaptic branches distributed according to the conduction delay
float	*s[NMAX+2][D];				// Synaptic strength
float	*sd[NMAX+2][D];				// Synaptic derivative
float	*sds[NMAX+2][D];			// ST-STDP strength
float	DA = 0.0;						// Extracellular dopamine
float DA_decay_rate = 0.0;
float LR = 0.0;						//Learning Rate, RGM
float sd_decay_rate = 0.951;//default decay rate, down to zero at ~6 seconds

#ifdef FindGlobalAxons
short	S_local[Br]; 				// synapses 0-S_local are local, synapses S_local -> N_pre are global
#endif


const	float	E_glu	=	0.0;			// Reverse potential for NMDA and AMPA currents
const	float	E_GABAa	=	-70.0;
const	float	E_GABAb	=	-90.0;
const	float	E_GABAc	=	-90.0;

int		Ntypes=0;					// The number of different neuronal types; to be read from file.



// The percentage of different neuronal types
float	Types[MAX_TYPES] = {0.0};  // set them all to zero; types are read in from connect.dat 

int			itype[MAX_AREAS][MAX_TYPES+1];
char		AreaNames[MAX_AREAS][128];
int			type_count[MAX_AREAS][MAX_TYPES];  // Number of neurons of each type.
group		group_data[MAX_AREAS][MAX_TYPES];	// Everything you need to know about a particular type of neuron in a particular area.

vector<string> layerNames;

// ********************************
// The following data structures define the repertoire of neurons.  These params will be read from the neuronspecs.dat
//
vector<string> Names;
vector<int> Neuronal_Class;
vector<int> glutamatergic;
vector<int>	soma_layer;	// 1 through 6 for layers 1 through 6 in cortex + 7-thalamus,
vector<float>	C;
vector<float>	K;
vector<float>	vr;
vector<float>	vt;
vector<float>	spike_peak_soma;
vector<float>	spike_peak_dend;	
// Up: from daughters, Down: from mother
//DEBUG
vector<float>	Gdendrites_Up ; 	
vector<float>	Gdendrites_Down ;
vector<float>	a;
vector<float>	b;
vector<float>	spike_reset_c_soma;	
vector<float>	spike_reset_c_dend;	
vector<float>	spike_reset_d_soma;	
vector<float>	spike_reset_d_dend;	
vector<float>	U_upper_bound;

vector<float> Background_AMPA_for_type;
vector<float> Background_GABA_for_type;
vector<float> decay_ampa;
vector<float> decay_nmda;
vector<float> decay_gabaa;
vector<float> decay_gabab;

// Axonal RADIUS in mm:  NOT USED
vector<vector<float> >	axons; 

// added 6/1/11 AEK
vector<double> gabac_const;
vector<double> decay_gabac;
vector<int> gabac_start_msec;
vector<int> gabac_stop_msec;

//*******************************************

float	frate[MAX_TYPES]; 

// distribution of synapses for different layers for different neuronal types
int	synapses[MAX_AREAS][MAX_TYPES][Layers]={}; // Read in from microcircuitry.dat

vector<vector<connection> >	SYN[MAX_AREAS][MAX_TYPES][Layers];  // To save space, all vectors of synaptic contact percentages are zero, unless microcircuitry.dat has
																// connections specified in that area.
/*
 Anatomy matrix. Cortico-cortical connections, brainstem, and sensory connections are in the last column.  The last two are for future use.
 The file that fills this structure takes on the following format, similar to the table in Izhikevich and Edelman, 2008 PNAS paper, supporting information.
 
 posttype layer % cells synapses (%synapses from each presynaptic neuron type)
 nb1	L1		1.5		8890	10.1,  6.3,  0.6,  1.1,  0.0,  0.0,  0.1,  0.0,  0.0,  0.1,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  4.1,  0.0,  0.0,  0.0, 77.6,  0.0,  0.0, 
 p23	L23		26.0	5800	 0.0, 59.9,  9.1,  4.4,  0.6,  6.9,  7.7,  0.0,  0.8,  7.4,  0.0,  0.0,  0.0,  2.3,  0.0,  0.0,  0.8,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
 p23	L1		0		1306	10.2,  6.3,  0.1,  1.1,  0.0,  0.0,  0.1,  0.0,  0.0,  0.1,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  4.1,  0.0,  0.0,  0.0, 78.0,  0.0,  0.0,
 ...
 */




// Neuronal Classes: 
// other=0, (ss,p)=1, b=2, nb=3, TC=4, T(inh)=5, 
const short Neuronal_classes = 6;

// Eugene's original: const	float	sm[Neuronal_classes]={6.0, 10.0, 6.0, 4.0, 20.0, 5.0};	
// sm[pre_class][post_class]


// maximal synaptic weight/conductance on different types of neurons 
// nb1=14, p23=250


// Short-term synaptic plastity (STSP) types: 
// 0 - no plasticity,
// 1 RS->RS, RS->FS, FS->RS, FS->FS (i.e., p->p, p->b, b->p, b->b)
// 2 RS->LTS (i.e., p->nb)
// 3 TC->RS  (i.e., tc->p)
// 4 TC->FS	 (i.e., tc->b)

const int	STSP_types=5;

//float	STSP_p[STSP_types]=		{1, 	0.6,	1.5,	0.7,	0.5}; 	
float	STSP_p[STSP_types]=		{1, 	0.8,	1.5,	0.7,	0.5}; 	
float STSP_tau[STSP_types]=	{1,	150,	100,	150,	200};

const int	ST_plast=1000;				// Length of STSP window of interest
float	STSP[STSP_types][ST_plast];	// Pre-computed exponents to determine the amount of recovery

short STSP_Type[Neuronal_classes][Neuronal_classes] = {
	// other		sp,p		b		nb		TC		T(inh)
	0,				0,			0,		0,		0,		0,			// other->
	0,				1,			1,		2,		0,		0,			// p-> (i.e., RS->)
	0,				1,			1,		0,		0,		0,			// b-> (i.e., FS->)
	0,				0,			0,		0,		0,		0,			// nb-> (i.e., LTS->)
	0,				3,			4,		0,		0,		0,			// TC->
	0,				0,			0,		0,		0,		0};		// T(inh)->


short STSP_s_var[STSP_types] = {0, 0, 1, 0, 1}; // the identity of s variable for each type of STSP_s
short STSP_s_type[2][Neuronal_classes] = {  // i.e., which p and tau to use for each STSP_s
	0, 1, 1, 0,  3, 0,	//  first variable in STSP_s
	0, 2, 0, 0,  4, 0};	// second variable in STSP_s

// How this works while EPSPs and IPSPs are delivered:
// neuron number -> neuronal class;
// pre neuronal class + post neuronal class -> STSP type
// STSP type - > s variable

// How this works when a pre neuron fires
// neuron number -> neuronal class
// neuronal class -> the type of plasticity of s variables (via STSP_s_type)
// update STSP_s according to the appropriate STSP_p and STSP_tau

//int	training_pts;;
int APEX_button_press = 0;
int APEX_REWARD = 0;

int	read_training_array( char* fname, const int training_dim )
{
	
	int i=0;
	int k=0;
	
	//float array[100000][2];
	float x;
	//char* fname="force_dv.dat";
	fstream gfile(fname);
	while(gfile>>x)
	{
		if (i==0) {
			training_array[k][0]=x;
			i=1;
		}
		else if (i==1) {
			training_array[k][1]=x;
			i=2;
		}
		else {
			training_array[k][2]=x;
			i=0;
			k++;
		}
	}
	
	//training_pts=k;
	return k;
	//cout << "finishe input training data. " << k/training_dim << "\n";
	//cout << "finish input training data. " << k << " " << training_pts << "\n";
}


// Initialize Branch
void	Make_Branch(int i_neu, unsigned char L, int Down, int Up, int Sister)
{
	
	Br_N[N_br]=i_neu;
	Br_L[N_br]=L;
	
	Br_Down[N_br]=Down;				// The link to the mother branch
	Br_Up[N_br]=Up;					// The link to the daughter branch
	Br_Sister[N_br]=Sister;			// The link to the sister branch
	
	
	if (N_br >= Br-1) 
	{
		fprintf(stderr,"Allocation error on proc=%d: Neuron %d needs more branches. Allocated: %d. Max: %d. Increase 'Br' \n", myid, i_neu, N_br, Br);
		exit(0); 
	}
}


/*
 Neuronal dendritic tree (move from botton to up)
 
 789	next layer					9 - apical
 \ 
 456	next layer					6 - apical
 \ 
 23		2 - basal dendrite,	3 - apical
 \ 
 s		soma
 
 */

void	Make_Neuron(int i_neu, short type, short area) 
{	// creates dendritic structure for the neuron
	
	unsigned char	layer_id;
	int		Mother, sister;
	int		ignore_soma;
	
	for (layer_id=0;layer_id<Layers; layer_id++)
		if (synapses[area][type][layer_id]>0) // skip if no synapses are received in this layer
		{
			ignore_soma=0;
			if (Soma[i_neu]<0) 	  // still need to allocate the soma compartment
								  // (assumes that there are no synapses below the soma)
			{
				Soma[i_neu]=N_br; // Somatic branch is like any other branch of the tree
				Mother=N_br;
				Make_Branch(i_neu, layer_id, -1, -1, -1);
				N_br++;
				ignore_soma=1;
				if (ceil(synapses[area][type][layer_id]/M) > 1) 
					Br_Up[Mother]=N_br; // there will be more branches  
			}
			else
				Br_Up[Mother]=N_br; // there will be more branches  
			
			
			
			
			// Now, add branches to this layer
			for (sister=ignore_soma; sister<ceil(synapses[area][type][layer_id]/M); sister++)
			{
				Make_Branch(i_neu, layer_id, Mother, -1, N_br+1); // this branch=N_br, next one=N_br+1
				N_br++;
			}
			// Now, N_br is looking at the next (empty) branch
			Mother=N_br-1; // the last branch is the apical part, from which the rest of the tree grows
			Br_Sister[Mother]=-1; // no more syster branches in this layer
		}
	
	if (Soma[i_neu]<0) 	  
	{
		// still need to allocate the soma compartment.  This means the cell has no synapses.  
		// Allocate the neuron anyway with one compartment to allow for dummy neurons.
		for (layer_id=0;layer_id<Layers; layer_id++){
			if(SYN[area][type][layer_id].size() > 0){
				// The first layer mentioned.  This is where the cell body should go.
				
				Soma[i_neu]=N_br; // Somatic branch is like any other branch of the tree
				Mother=N_br;
				Make_Branch(i_neu, layer_id, -1, -1, -1);
				N_br++;
			}
		}
	}
	
}




void	allocate_branches()
{
	int i;
	
	
	for (i=myid;i<N;i=i+NP) 
	{ // for each neuron in this subnetwork
		if(( myid == watchid) && (i%1000<NP)) fprintf(stderr,"Type %d, neuron %d on processor %s\n", Type[i], i, processor_name);
		Make_Neuron(i, Type[i], Area[i]);
	}
	
	if( myid == watchid) fprintf(stderr,"\nUsed compartments/branches: %d \n", N_br);
	
}


int	MoveAlongApicalDendrite(int br)
{
	if (Br_Up[br] < 0)
		return -1; // this is the last branch.
	
	br = Br_Up[br]; // one level up
	
	while (Br_Sister[br] > 0)   // i.e., there is another one.
		br = Br_Sister[br];		// move right toward the apical branch, which is the last
	
	return br; // the last sister is a segment of the apical dendrite.
	
}

// Returns -1 when no more compartments for this neuron, otherwise returns next br number.
// Assumes you start at soma and loop through comparments:
// br=Soma[nrn];
// while(br>=0){
//	... 
//	br=NextCompartment(br);
// }
int NextCompartment(int br){
	if(br < 0) return -1;  // This indicates that the Soma[nrn] is not onboard.  Obviously no more compartments.
	if(Br_Sister[br]>=0) return Br_Sister[br];
	// If no more sisters, check to see if we have daughters
	if(Br_Up[br]>=0) 
		return Br_Up[br];
	else
		// We're' done.
		return -1;
}



areas area_data;
int N_areas;

//**************************************************************
// init AreaNames and N_areas and area_data.
void	read_areas(){
	
	// Read the area names and 
	
	area_data.load_areas("areas.dat");
	
	for(int i = 0; i < area_data.get_N_areas();i++){
		string area_name = area_data.get_area_name(i);
		strcpy( &(AreaNames[i][0]), area_name.c_str());
		AreaNames[i][area_name.size()]='\0';  // EOS character.
	}
	N_areas = area_data.get_N_areas();
	
	
}


// returns density in neurons/mm^2 of a projection from pre_area,pre_type to post_area.
// 
inline float projection_density(int post_area,int pre_area,int pre_type){
	
	float pre_neurons = group_data[pre_area][pre_type].rows*group_data[pre_area][pre_type].cols;
	int anytype = 0; // We don't care which type;  all post_area neuron types have the same area size (stored redundantly for convenience.)
	float post_mm_squared = group_data[post_area][anytype].area_width_mm* group_data[post_area][anytype].area_width_mm;
	
	return pre_neurons/post_mm_squared;
	
}




// Used to understand the network;
// Pre: call after allocate_synapeses();
// Post: file "connection_prob.dat" is written;

void write_connection_probabilities(string area_name){
	
	string s = "connection_prob_" + area_name + ".dat";
	
	ofstream fout(s.c_str());
	
	int post_area = area_data[area_name].area_number;
	fout << area_name << ":"<< post_area << endl;
	
	
	//	for (int pre_type=0;pre_type<Ntypes;pre_type++){ // the last element, Ntypes,Ntypes+1,Ntypes+2 correspond to cortico-cortical, brainstem, and sensory synapses, respectively
	//		fout << setw(8) << Names[pre_type] << " ";
	//	}
	fout << endl;
	fout << "post:               ";
	for (int post_type = 0; post_type<Ntypes; post_type++){
		
		for( int layer = 0; layer < Layers; layer ++ ){
			if( SYN[post_area][post_type][layer].size() > 0 ){  
				// synapses in that layer.
				fout << setw(8) << Names[post_type] << " " << setw(4) << layerNames[layer] << " " << setw(6) << type_count[post_area][post_type];
				for (int pre_area=0;pre_area<N_areas;pre_area++){
					for (int pre_type=0;pre_type<Ntypes;pre_type++) 
					{
						float p_synapse;
						int Nsyn = synapses[post_area][post_type][layer]*0.01*SYN[post_area][post_type][layer][pre_area][pre_type].pct_pre_synapses;	// the number of synapses needed
						if( Nsyn > 0 ) {
							float dmax2 = SYN[post_area][post_type][layer][pre_area][pre_type].max_radius_mm/Rc;		// normalization  - dmax2 is radius in mm.
							dmax2 = dmax2*dmax2; // squared
							
							float density = projection_density(post_area,pre_area,pre_type);
							float neurons_in_reach = density * pi*dmax2;
							if( neurons_in_reach < 0.1) {
								p_synapse = 0.0;
							}
							else{
								// Don't divide by zero. 
								p_synapse = Nsyn/neurons_in_reach;
							}
							fout << "from: " << AreaNames[pre_area] << Names[pre_type] << " " << setw(8) << fixed << setprecision(5) << p_synapse << "; ";
						}
					}
				}
				fout << endl;
			}
		}
	}
	fout.close();
	
}




//***************************************************************
void	allocate_somas_on_cortex()	   
{
	int	i, j, rj, t, type, cell_type;
	float x1, x2, x3;
	
	
	
	neuronspecs ns;
	ns.read_neuronspecs();
	Ntypes = ns.get_num_types();
	layerNames = ns.get_layer_names(); 
	
	ns.get_all(       
			   Names,
			   Neuronal_Class,
			   glutamatergic,
			   soma_layer,	
			   C,
			   K,
			   vr,
			   vt,
			   spike_peak_soma,
			   spike_peak_dend,	
			   Gdendrites_Up , 	
			   Gdendrites_Down ,
			   a,
			   b,
			   spike_reset_c_soma,	
			   spike_reset_c_dend,	
			   spike_reset_d_soma,	
			   spike_reset_d_dend,	
			   U_upper_bound,
			   Background_AMPA_for_type,
			   Background_GABA_for_type,
			   decay_ampa,
			   decay_nmda,
			   decay_gabaa,
			   decay_gabab,
			   axons,
			   gabac_const,
			   decay_gabac,
			   gabac_start_msec,
			   gabac_stop_msec
			   );
	
	if( myid==watchid){
		printvector(Names);
		printvector(Neuronal_Class);
		printvector(glutamatergic);
		printvector(soma_layer);	
		printvector(C);
		printvector(K);
		printvector(vr);
		printvector(vt);
		printvector(spike_peak_soma);
		printvector(spike_peak_dend);	
		printvector(Gdendrites_Up ); 	
		printvector(Gdendrites_Down );
		printvector(a);
		printvector(b);
		printvector(spike_reset_c_soma);	
		printvector(spike_reset_c_dend);	
		printvector(spike_reset_d_soma);	
		printvector(spike_reset_d_dend);	
		printvector(U_upper_bound);
		printvector(Background_AMPA_for_type);
		printvector(Background_GABA_for_type);
		printvector(decay_ampa);
		printvector(decay_nmda);
		printvector(decay_gabaa);
		printvector(decay_gabab);
		// added 6/1/11 AEK
		printvector(gabac_const);
		printvector(decay_gabac);
		printvector(gabac_start_msec);
		printvector(gabac_stop_msec);
	}		
	
	
	// Create connect.dat and groups.dat before we start; info is read from areas.dat and microcircuitry.dat 
	// written by matlab script make_microcircuitry.m in the matlab subdirectory.
	neurons neuron_data;
	bool write_flag = (myid == watchid);
	neuron_data.set_positions(group_data, x, write_flag);
	
	read_areas();
	
	
	microcircuitry<MAX_AREAS,MAX_TYPES,Layers> micro;
	micro.read_microcircuitry(myid);
	if(myid==watchid) cerr << "*** after read_microcircuitry *** " << ceil(synapses[0][0][0]/M) << endl;
	//if(myid == watchid) micro.print();
	
	micro.get_microcircuitry(synapses, SYN);
	
	if(myid==watchid) cerr << "*** after get_microcircuitry *** " << ceil(synapses[0][0][0]/M) << endl;
	if(myid == watchid) {
		ofstream fout("microcircuitry_used.dat");
		micro.print(fout,SYN);
		fout.close();
		
	}	
	
	int neuron_index=0;
	for(int area = 0; area < N_areas;area++){
		for(type=0;type<Ntypes;type++){
			int n = group_data[area][type].rows * group_data[area][type].cols;
			type_count[area][type] = n;
			
			Types[type]+=n;
			for( int i=0;i<n;i++){
				
				Type[neuron_index]=type;
				Area[neuron_index]=area;
				Soma[neuron_index]=-1;  // May not have a soma if neuron isn't on this processor.
				neuron_index++;
			}
		}
	}	
	
	
	if( neuron_index > NMAX ){		
		if( myid == watchid) fprintf(stderr,"Too many neurons (lines) in connect.dat.  MAX set to %d.\n", NMAX);
		exit(0);
	}
	N=neuron_index;
	
	cout << "Total Neurons: " << N << endl;
	
	// Set up the itype structure.
	int count = 0;
	for(int area = 0; area < N_areas; area++){
		itype[area][0]=count;  // pick up where we left off.
		for (j=1;j<Ntypes;j++)
		{
			itype[area][j] = itype[area][j-1]+type_count[area][j-1];
			count += type_count[area][j-1];
			cout << itype[area][j] << " " << endl;
		}
	}
	
	
	for (j=0;j<Ntypes;j++)
	{
		Types[j] = Types[j]/N*100.0;
	}
	
	//save data	for DEBUG
	if( myid == 0) {
		FILE * fs = fopen("neurons_debug.dat","w");
		for (i=0;i<N;i++)
			fprintf(fs, "%4.4f  %4.4f %4.4f \n", x[i][0], x[i][1], x[i][2]);
		fclose(fs);
	}
	
}



// Need this function when the number of required synapses is less than 1.
int	HowManySynapses( float x, int max_number)
{
	int ans;
	ans = (int)floor(x);
	if (random01() < x-ans)
		ans++;
	
	if (ans>max_number)
		ans = max_number;
	
	return ans;
}

// Now we need to convert pre to post, taking into account delay information in pre_delay.
// That is, go through all neurons and combine synapses  
// N+1 - brainstem, N+2 - sensory

static 	int		i_delays[NMAX+2][D];			// The current index in the distribution of conduction delays from pre neuron to all post neurons
void Convert_pre_to_post()
{
	int	i,j, br, p, grand_total;
	short	del;
	grand_total = 0;
	
	if( myid == watchid) fprintf(stderr,"Converting pre to post...\n");
	
	for (i=0;i<N+2;i++) 
		for (del=0; del<D; del++)
		{
			delays[i][del]=0;
		}
	
	if( myid == watchid) fprintf(stderr,"Counting synapses...\n");
	
	for (br=0;br<N_br;br++)
		for (j=0;j<M;j++)
			if (pre[br][j]>=0)  // a valid index
			{
				p = pre[br][j];
				del = pre_delay[br][j];
				delays[p][del]++;
				grand_total++;
			}
	
	if( myid == watchid) fprintf(stderr,"Total synapses= %d\n",grand_total);
	
	
	if( myid == watchid) fprintf(stderr,"Allocating memory...\n");
	for (i=0;i<N+2;i++)
	{
		//		if ( myid == watchid) 
		if (i%10000==0) fprintf(stderr,"myid=%d, Neuron %d\n",myid, i);
		
		for (del=0; del<D; del++)
		{ 
			i_delays[i][del]=0;
			post[i][del]=	(int *)	 malloc(delays[i][del]*sizeof(int));
			s[i][del]=		(float *)malloc(delays[i][del]*sizeof(float));
			sd[i][del]=		(float *)malloc(delays[i][del]*sizeof(float));
			sds[i][del]=	(float *)malloc(delays[i][del]*sizeof(float));
			
			if ((post[i][del]==NULL) | (s[i][del]==NULL) | (sd[i][del]==NULL) | (sds[i][del] == NULL))
			{
				fprintf(stderr,"Cannot allocate memory for neuron %d\n, processor myid=%d",i, myid);
				exit(0);
			}
		}
	}
	if( myid == watchid) fprintf(stderr,"Making Connections...\n");
	for (br=0;br<N_br;br++)
		for (j=0;j<M;j++)
			if (pre[br][j]>=0) // a valid index
			{
				p = pre[br][j];	// pre-synaptic neuron
				del = pre_delay[br][j];  
				i=i_delays[p][del]++;
				post[p][del][i]=br;
				s_pre[br][j]=&s[p][del][i];
				sd_pre[br][j]=&sd[p][del][i];
				sds_pre[br][j]=&sds[p][del][i];
				
			}
	
	// check point: i_delays should be equal to delays at the end
	
	if( myid == watchid) fprintf(stderr,"Confirming initialization.\n");
	for (i=0;i<N+2;i++)
		for (del=0; del<D; del++)
			if ((delays[i][del] !=i_delays[i][del]))
			{ 
				//neuron 650764, delay 16, myid=53
				fprintf(stderr,"Initialization error: neuron %d, delay %d, myid=%d",i,del, myid);
				exit(0); 	
			}
	
}

// find the delay from presynaptic neuron neu to the postsynaptic branch br. 
// Delay is in [0, D-1] range
// Cortico-cortical Connection could be 
//		global (myelinated) maximal: D/2, corresponding to the maximal axonal length of 2Rc
//		local (nonmyelinated). Maximal: D/2	corresponding to longest axonal span of 1.2 mm
//		layer-to-layer = 1 ms for each layer difference (corresponds to .1 meter/second or .1mm/msec (Salin&Prince, 1996). 
//		delays from thalamus to cortex = 1 ms + local layer-to-layer 
//		delays from cortex to thalamus
//			layer 5 to cortex = 1 ms 
//			lyaer 6 to cortex = random [0,D-1] per neuron

const	float	unmyelinated_speed = 1.2/(D/2);

short	Determine_Delay(short type_pre, short type_post, float d, float long_axon) // the distances d,long_axon are squared (rescaled distances on the unit sphere)
{
	int	p;
	short	layer_pre, layer_post;
	short del;
	
	layer_pre = soma_layer[type_pre];
	layer_post= soma_layer[type_post];
	
	if (layer_post==7) // thalamic neuron
	{
		if (layer_pre==5)
			return 1 + (short) ceil(sqrt(d)*Rt/unmyelinated_speed); // 1 ms + axonal span in the thalamus
		return getrandom(D);	// error: different delays for the same layer 6 pre neuron
	}
	
	if (layer_pre==7) // input from thalamus
	{
		del= 1+1*(6-layer_post)+ (short) ceil(sqrt(d)*Rc/unmyelinated_speed);	// 1 ms + 1ms for each cortical layer + axonal span in cortex
		
		if (del >= D)
			del = D-1;
		
		return del;
	}
	
	// Both are cortical neurons
	
	del =  abs(layer_pre-layer_post); // layer-to-layer delay
	del += (short) floor( sqrt(d)*Rc/unmyelinated_speed + long_axon*Rc/myelinated_speed); // local + global delays
	
	if (del >= D)
		del = D-1;
	
	
	return del;
}

//********************************
// Find the area under the normal curve with sigma, sigma_mm, from 0 to x.
// Approximation to erf (Wikipedia). 
float norm_integral(float sigma_mm, float x){
	
	x= x/(sigma_mm*sqrt(2)); // adjust for normal distribution.
	
	// approximation to erf(x)
	const float a=-8*(pi-3)/(3*pi*(pi-4));
	
	return sqrt(1-exp(-x*x*(4/pi+a*x*x)/(1+a*x*x) ));
	
}



//********************************
float gaussian_pdf(float sigma, int Nsyn, int neurons_in_reach, float d_mm, float max_d_mm ){
	
	bool uniform;
	float g;
	
	float scale = (float)Nsyn*(2*max_d_mm)/norm_integral(sigma,max_d_mm)/neurons_in_reach;
	
	uniform = false;
	if(scale > sigma*sqrt(2*pi))
	{
#ifdef WARNINGS
		if(myid == watchid)
			cerr << "warning: not enough cells in reach to get gaussian pdf; switching to uniform distribution" << endl;
#endif
		uniform = true;
	}
	
	if (uniform){
		float one_over_dmax = 1.0/max_d_mm;
		g= (float)Nsyn/neurons_in_reach*(1.6667-d_mm*one_over_dmax);
	}
	else{
		float a = 1.0/(sigma*sqrt(2*pi));
		g= scale*a*exp(-(d_mm*d_mm)/(2*sigma*sigma));
	}
	
	
	return g;
	
}

// If a cell in area A projects to a cell in area B, we need to find the topographic projection point, <axonx,axony> and axon length, dist.
void translate_pre_to_post_coord(int pre_id, int post_id, float &axonx, float &axony, float &dist){
	
	int anytype = 0; // Don't care which type since they all have the same area info.
	
	int post_area = Area[post_id];
	int pre_area  = Area[pre_id];
	
	float post_area_width_mm = group_data[post_area][anytype].area_width_mm;
	float pre_area_width_mm  = group_data[pre_area ][anytype].area_width_mm;
	float scale = post_area_width_mm/pre_area_width_mm;
	
	float x = ::x[pre_id][0];
	float y = ::x[pre_id][1];
	
	
	float tx0 = group_data[post_area][anytype].x0;
	float ty0 = group_data[post_area][anytype].y0;
	
	float x0 = group_data[pre_area][anytype].x0;
	float y0 = group_data[pre_area][anytype].y0;
	
	// Rescale the relative position of this neuron from the area origin (x0,y0) and add to the target area origin (tx0,ty0)
	axonx=tx0 + (x-x0)*scale;
	axony=ty0 + (y-y0)*scale;
	
	// Now that we know the coordinates of this pre_neuron's axon in the target area, calculate the distance of this white-matter axon.
	
	float tx = ::x[post_id][0];
	float ty = ::x[post_id][1];
	
	float dx = (axonx-tx);
	float dy = (axony-ty);
	dist= sqrt( dx*dx + dy*dy );
	
}


void allocate_synapses()
{
	int i,j, k, br, post_neuron, start_search, candidate;
	short post_type, pre_type;
	unsigned char layer;
	short Nsyn;
	float d2, dmax2, R, d_long;
	int	p;
	
	
	
	for (br=0;br<N_br;br++)
	{
		
		post_neuron=Br_N[br];
		post_type=Type[post_neuron];
		int post_area = Area[post_neuron];
		
		// assert:  SYN[post_area][post_type][layer].size() > 0 here since we are looking only at branches that really exist, and all branches live in a post_area,
		//    have a post_type and exist in a given layer.
		
		float post_area_width_mm=area_data[ AreaNames[post_area] ].width_mm;
		
		layer=Br_L[br];
		N_pre[br]=0;	// the number of allocated synapses
		if (layer < 5 ) 
			R=Rc;	// we are in cortex
		else		
			R=cortex_to_thalamus_size_ratio;	// we are in thalamus
		
		
		for (int pre_area=0;pre_area<N_areas;pre_area++){
			
			for (pre_type=0;pre_type<Ntypes;pre_type++) // the last element, Ntypes,Ntypes+1,Ntypes+2 correspond to cortico-cortical, brainstem, and sensory synapses, respectively
			{
				
				// number_wanted is the number of synapses for this cell type in this compartment;  cases are different if M is very large to make things psuedo-single compartment or if M << synapses, which is the normal multi-compartment model.  That is, if M is huge then we'll accidentally allocate way too many synapses under the old method
				// JLM JGF 6/2/10
				
				int number_wanted = (int)min((double)M, ceil( synapses[post_area][post_type][layer]) );
				Nsyn = HowManySynapses(0.01 * number_wanted * SYN[post_area][post_type][layer][pre_area][pre_type].pct_pre_synapses, number_wanted-N_pre[br]);	// the number of synapses needed
				
				if( Nsyn < 1 ) continue;
				
				// We know there are synapses between these groups.
				float pre_area_width_mm=area_data[ AreaNames[pre_area] ].width_mm;
				
				connection conn = SYN[post_area][post_type][layer][pre_area][pre_type];
				
				float dmax = conn.max_radius_mm*R;		// normalization  - dmax2 is radius in mm.
				dmax2 = dmax*dmax; // squared, because d2 below is squared
				
				float dmin = conn.min_radius_mm*R;		// normalization  - dmax2 is radius in mm.
				float dmin2 = dmin*dmin; // squared, because d2 below is squared
				float sigma = conn.sigma_mm*R;		// don't need to normalize if sigma and dmax are unique to pre/post?
				
				
				//Instead of determining the number of synapses; determine the prob of connection with
				// this cell type.  Loop over all neurons looking for synapses.
				
				float density = type_count[pre_area][pre_type]/(post_area_width_mm*post_area_width_mm);
				
				float neurons_in_reach = density * (pi*dmax2 - pi*dmin2);  // Assumes dmin2 will be 0 for gaussian, and nonzero for surround projections.
				
				if( neurons_in_reach<0.1)
				{
#ifdef WARNINGS
					if(myid==watchid)
						cerr << "Warning: no neurons in reach for type: " << Names[post_type] << " from type " << Names[pre_type] << "!\n";
#endif
					continue; // Don't divide by zero. 
				}
				float p_synapse = Nsyn/neurons_in_reach;

#ifdef WARNINGS				
				if( (myid == watchid) && (p_synapse>1) )
					cerr << "Warning: not enough neurons in reach for type: " << Names[post_type] << " from type " << Names[pre_type] << "in layer" << layer << "area" << post_area <<". neurons_in_reach: " << neurons_in_reach << " " << "pi*dmax^2: " << pi*dmax2<< "; density = " << density << "p_synapse = " << p_synapse << endl;
#endif		
		
				//if( myid == watchid && (br%100 == 0)) fprintf(stderr,"P(synapse)=%f for neuronal type %d from neuronal type %d \n{0=nb1, 1=p23, 2=b23, 3=nb23, 4=ss4L4, 5=ss4L23, 6=p4, 7=b4, 8=nb4, 9=p5L23, 10=p5L56, 11=b5, 12=nb5, 13=p6L4, 14=p6L56, 15=b6, 16=nb6, 17=TCs, 18=TCn, 19=TIs, 20=TIn, 21=TRN}\n in layer type %d\n {0=L6, 1=L5, 2=L4, 3=L23, 4=L1, 5=Thal, 6=RTN} \n",p_synapse, type, pre_type, layer);
				// We need: 
				// Nsyn synapses 
				// from pre_type neurons 
				// within sqrt(dmax2) mm of this neuron (located at x[neuron])
				
				
				k=0;
				for(candidate=itype[pre_area][pre_type];candidate<itype[pre_area][pre_type+1]&&N_pre[br]<M;candidate++){
					
					float tx,ty;
					float tz=0.0f;
					if( pre_area!=post_area ){
						// Global connection.  Translate into local coordinates.
						translate_pre_to_post_coord(candidate, post_neuron, tx, ty, d_long);
					}
					else{
						// Areas are the same; no coordinate conversion necessary, and no between-area delays.
						tx=x[candidate][0];
						ty=x[candidate][1];
						d_long = 0.0f;
					}
					
#ifdef WRAPAROUND
					d2 = (min(fabs(x[post_neuron][0]-tx),post_area_width_mm-fabs(x[post_neuron][0]-tx))*min(fabs(x[post_neuron][0]-tx),post_area_width_mm-fabs(x[post_neuron][0]-tx)))
					+(min(fabs(x[post_neuron][1]- ty),post_area_width_mm-fabs(x[post_neuron][1]-ty))*min(fabs(x[post_neuron][1]-ty),post_area_width_mm-fabs(x[post_neuron][1]-ty)))
					+(min(fabs(x[post_neuron][2]- tz),post_area_width_mm-fabs(x[post_neuron][2]-tz))*min(fabs(x[post_neuron][2]-tz),post_area_width_mm-fabs(x[post_neuron][2]-tz)));
#else
					d2 = (x[post_neuron][0]-tx)*(x[post_neuron][0]-tx) + (x[post_neuron][1]-ty)*(x[post_neuron][1]-ty) + (x[post_neuron][2]-tz)*(x[post_neuron][2]-tz);
#endif
					
					//if( d2 < dmax2 && random01() < p_synapse*(1.6667-sqrt(d2*one_over_dmax2)) ){ // higher p(conn) if closer.
					if( d2 < dmax2 ){
						if( dmin>.0001 ){
							if( d2 >= dmin2 ){
								// Surround projection.  Transform it into a Gaussian shifted s.t. the center of the Gaussian is midway between dmin and dmax.
								float center_radius_mm = (dmin+dmax)*0.5f;
								float dist_from_center = center_radius_mm - sqrt(d2);
								float dmax_surround = center_radius_mm-dmin;
								
								if( random01() < gaussian_pdf( sigma, Nsyn, neurons_in_reach, dist_from_center, dmax_surround  ) ){ // higher p(conn) if closer.
									pre_delay[br][N_pre[br]]=Determine_Delay(pre_type, post_type, dist_from_center*dist_from_center, d_long); // square the distance since function takes sqrt.
									pre[br][N_pre[br]++]=candidate;		
									k++;
								}
							}
						}
						else{
							// Standard Gaussian projection.
							if( random01() < gaussian_pdf( sigma, Nsyn, neurons_in_reach, sqrt(d2), dmax  ) ){ // higher p(conn) if closer.
								pre_delay[br][N_pre[br]]=Determine_Delay(pre_type, post_type, d2, d_long); // TODO: assign delays to neurons, not branches
								pre[br][N_pre[br]++]=candidate;		
								k++;
							}
						}
					}
				}
				//d2 = (x[neuron][0]-0.5)*(x[neuron][0]-0.5) + (x[neuron][1]-1.5)*(x[neuron][1]-1.5) ;
				//if(myid==watchid && d2<.5*.5 )cout << Nsyn << " " << k << " " << pre_type << " " << N_pre[br]<< endl;
			}
		}
		if (( myid == watchid) && (br%1000==0)) fprintf(stderr,"Synapses for branch: %d. Neuron:%d, layer:%d,=%d \n", br, post_neuron, layer, N_pre[br]);
		
		
#ifdef FindGlobalAxons
		S_local[br]=N_pre[br];
#endif
		
		// NOTE: no special case for cortico-cortical connections any longer to allow for arbitrary type connections between groups to allow feedforward
		//  and feedback and lateral type reentry between areas.  Also, this allows individual cells to project to multiple cortical areas.  Finally, 
		//  this allows independent normalization of different afferent sources.
		
		// Brainstem
		/*
		 pre_type=Ntypes+1;
		 {
		 Nsyn = HowManySynapses(0.01*M*SYN[row][pre_type], M-N_pre[br]);	// the number of synapses needed
		 for (k=0;k<Nsyn;k++)
		 {
		 pre_delay[br][N_pre[br]]=0;
		 pre[br][N_pre[br]++]=N; // Brainstem input;	
		 }
		 }
		 
		 // Sensory
		 pre_type=Ntypes+2;
		 {
		 Nsyn = HowManySynapses(0.01*M*SYN[row][pre_type], M-N_pre[br]);	// the number of synapses needed
		 for (k=0;k<Nsyn;k++)
		 {
		 pre_delay[br][N_pre[br]]=0;
		 pre[br][N_pre[br]++]=N+1; // sensory input; will be allocated later	
		 }
		 }
		 */
		
		for (k=N_pre[br];k<M;k++)
			pre[br][k]=-1; // in case some are unused
	}
	
    Convert_pre_to_post();
	
	
}

/*
 void save_allocated_somas(char fname[30])
 {
 FILE	*fs;
 fs = fopen(fname, "wb");
 fwrite(x,sizeof(float),3*N, fs);
 fwrite(Gpre,sizeof(int),N, fs);
 fwrite(&N_G,sizeof(int),1, fs);
 fclose(fs);
 }
 
 void load_allocated_somas(char fname[30])
 {
 FILE	*fs;
 fs = fopen(fname, "rb");
 fread(x,sizeof(float),3*N, fs);
 fread(Gx,sizeof(float),3*N, fs);
 fread(Gpre,sizeof(int),N, fs);
 fread(xydist,sizeof(short),N, fs);
 fread(&N_G,sizeof(int),1, fs);
 fclose(fs);
 }
 */

void save_allocated_branches(char fname[30])
{
	FILE	*fs;
	fs = fopen(fname, "wb");
	
	fwrite(&N_br,sizeof(int),1, fs);
	
	fwrite(Br_N,sizeof(int),Br, fs);
	fwrite(Br_L,sizeof(unsigned char),Br, fs);
	
	fwrite(Br_Down,sizeof(int),Br, fs);
	fwrite(Br_Up,sizeof(int),Br, fs);
	fwrite(Br_Sister,sizeof(int),Br, fs);
	
	fwrite(Soma,sizeof(int),N, fs);
	fwrite(Type,sizeof(unsigned short),N, fs);
	fwrite(Area,sizeof(unsigned short),N, fs);
	
	fwrite(Types,sizeof(float),MAX_TYPES, fs);
	fwrite(itype,sizeof(int),MAX_AREAS*(Ntypes+1), fs);
	
	fclose(fs);
}

void load_allocated_branches(char fname[30])
{
	
	FILE	*fs;
	int ret = 0;
	fs = fopen(fname, "rb");
	
	ret = fread(&N_br,sizeof(int),1, fs);
	
	ret = fread(Br_N,sizeof(int),Br, fs);
	ret = fread(Br_L,sizeof(unsigned char),Br, fs);
	
	ret = fread(Br_Down,sizeof(int),Br, fs);
	ret = fread(Br_Up,sizeof(int),Br, fs);
	ret = fread(Br_Sister,sizeof(int),Br, fs);
	
	ret = fread(Soma,sizeof(int),N, fs);
	ret = fread(Type,sizeof(unsigned short),N, fs);
	ret = fread(Area,sizeof(unsigned short),N, fs);
	
	ret = fread(Types,sizeof(float),MAX_TYPES, fs);
	ret = fread(itype,sizeof(int),MAX_AREAS*(Ntypes+1), fs);
	
	fclose(fs);
}

void save_allocated_synapses(char fname[30])
{
	FILE	*fs;
	fs = fopen(fname, "w");
	
	for(int i = 0; i<Br; i++){
		for(int j = 0; j < N_pre[i]; j++){
			if(pre[i][j] < N)fprintf(fs,"%d %d %f %d\n",Br_N[i],pre[i][j],*s_pre[i][j], pre_delay[i][j] );
		}
		//	fprintf(fs,"\n");
	}
	
	
	
	/*
	 fs = fopen(fname, "wb");
	 fwrite(N_pre, sizeof(short), Br, fs);
	 fwrite(pre,sizeof(int),Br*M, fs);
	 fwrite(pre_delay,sizeof(unsigned char),Br*M, fs);
	 */
	
	fclose(fs);
}


void save_sds(char fname[30])
{
	FILE	*fs;
	fs = fopen(fname, "w");
	
	for(int i = 0; i<Br; i++){
		for(int j = 0; j < N_pre[i]; j++){
			if(pre[i][j] < N)fprintf(fs,"%d %d %f \n",Br_N[i],pre[i][j],*sd_pre[i][j]);
		}
		//	fprintf(fs,"\n");
	}
	fclose(fs);
}	

void load_allocated_synapses(char fname[30])
{
	
	FILE	*fs;
	int ret = 0;	
	fs = fopen(fname, "rb");
	
	ret = fread(N_pre, sizeof(short), Br, fs);
	ret = fread(pre,sizeof(int), Br*M, fs);
	
	ret = fread(pre_delay,sizeof(unsigned char), Br*M, fs);
	
	Convert_pre_to_post();
	
	fclose(fs);
}

void save_dynamics(char fname[30])
{
	int i, del;
	FILE	*fs;
	int ret = 0;
	fs = fopen(fname, "wb");
	
	ret = fwrite(&sec,sizeof(int),1, fs);
	ret = fwrite(&B_firings,sizeof(int),1, fs);
	ret = fwrite(&N_firings,sizeof(int),1, fs);
	ret = fwrite(firings_time,sizeof(short),N_firings_max, fs);
	ret = fwrite(firings_index,sizeof(int),N_firings_max, fs);
	ret = fwrite(dendritic_spikes,sizeof(int),Br*sp_memory, fs);
	ret = fwrite(somatic_spikes,sizeof(int),N*sp_memory, fs);
	ret = fwrite(STSP_s, sizeof(float), N*2, fs);	
	ret = fwrite(Background_AMPA, sizeof(float), N, fs);	
	ret = fwrite(Background_GABA, sizeof(float), N, fs);	
	ret = fwrite(Br_AMPA,sizeof(float),Br, fs);
	ret = fwrite(Br_NMDA,sizeof(float),Br, fs);
	ret = fwrite(Br_NMDAVI,sizeof(float),Br, fs);
	ret = fwrite(Br_GABAa,sizeof(float),Br, fs);
	ret = fwrite(Br_GABAb,sizeof(float),Br, fs);
	ret = fwrite(Br_GABAc,sizeof(float),Br, fs);
	ret = fwrite(Br_V,sizeof(float),Br, fs);
	ret = fwrite(Br_U,sizeof(float),Br, fs);
	ret = fwrite(&DA,sizeof(float),1, fs); //
	ret = fwrite(&LR,sizeof(float),1, fs); //RGM: added learning rate
	
	for (i=0;i<N+2;i++)
		for (del=0;del<D;del++)
			fwrite(s[i][del], sizeof(float),delays[i][del], fs);
	
	for (i=0;i<N+2;i++)
		for (del=0;del<D;del++)
			fwrite(sd[i][del], sizeof(float),delays[i][del], fs);
	
	// write sds to file for ST-STDP
	for (i=0;i<N+2;i++)
		for (del=0;del<D;del++)
			fwrite(sds[i][del], sizeof(float),delays[i][del], fs);
	
	fclose(fs);
}

void load_dynamics(char fname[30])
{
	int i, del;
	
	FILE	*fs;
	int ret = 0;
	fs = fopen(fname, "rb");
	ret = fread(&start_sec,sizeof(int),1, fs);
	start_sec++; // We wrote it at the end of the loop, just before the counter was incremented.
	cout << "Start second." << start_sec << " **************************" << endl;
	ret = fread(&B_firings,sizeof(int),1, fs);
	ret = fread(&N_firings,sizeof(int),1, fs);
	ret = fread(firings_time,sizeof(short),N_firings_max, fs);
	ret = fread(firings_index,sizeof(int),N_firings_max, fs);
	ret = fread(dendritic_spikes,sizeof(int),Br*sp_memory, fs);
	ret = fread(somatic_spikes,sizeof(int),N*sp_memory, fs);
	ret = fread(STSP_s, sizeof(float), N*2, fs);
	ret = fread(Background_AMPA, sizeof(float), N, fs);	
	ret = fread(Background_GABA, sizeof(float), N, fs);	
	ret = fread(Br_AMPA,sizeof(float),Br, fs);
	ret = fread(Br_NMDA,sizeof(float),Br, fs);
	ret = fread(Br_NMDAVI,sizeof(float),Br, fs);
	ret = fread(Br_GABAa,sizeof(float),Br, fs);
	ret = fread(Br_GABAb,sizeof(float),Br, fs);
	ret = fread(Br_GABAc,sizeof(float),Br, fs);
	ret = fread(Br_V,sizeof(float),Br, fs);
	ret = fread(Br_U,sizeof(float),Br, fs);
	ret = fread(&DA,sizeof(float),1, fs);
	ret = fread(&LR,sizeof(float),1, fs); //RGM: added learning
	
	for (i=0;i<N+2;i++)
		for (del=0;del<D;del++)
			ret = fread(s[i][del], sizeof(float),delays[i][del], fs);
	
	for (i=0;i<N+2;i++)
		for (del=0;del<D;del++)
			ret = fread(sd[i][del], sizeof(float),delays[i][del], fs);
	
	// read sds from file for ST-STDP
	for (i=0;i<N+2;i++)
		for (del=0;del<D;del++)
			ret = fread(sds[i][del], sizeof(float),delays[i][del], fs);
	
	fclose(fs);
}




void dump_synapses(char fname[30])
{
	int br, j;
	FILE	*fs;
	
	if (myid != 0) return;
	
	fs = fopen(fname, "w");
	
	for (br=0;br<N_br;br++)
		if (glutamatergic[ Type[ Br_N[br] ] ] ==1)
			for (j=0; j<N_pre[br]; j++)
				if (glutamatergic[ Type[ pre[br][j] ] ] ==1)
					fprintf(fs, "%d %d %d %d %f3.3\n", Br_N[br], br, pre[br][j], Type[ pre[br][j] ], *s_pre[br][j]);
	fclose(fs);
}


//////////////////////////////////////////////////////////////////////////////
// If the weights become too high due to learning and too much depolarizing current
// is driving the cell, we should scale down the weight total on each pathway to avoid this.
// Pre: the max_depolarizing current for each cell onboard has been initialized
// Post: SYN[][][][][].init_max_sum has been adjusted if any post cell has too much current.
void homeostatic_scaling_adjustment(vector<float> &max_Isyn, bool glu_only=true){

	const float limit = -900;
	
	vector<vector<float> > maxmax(N_areas, vector<float>(MAX_TYPES,0.0f));
	
	// Find the max max_Isyn in each neuronal group 
	for(int post_nrn=myid;post_nrn<N;post_nrn+=NP){
		
		int post_type = Type[post_nrn];
		int post_area = Area[post_nrn];
		if( max_Isyn[post_nrn] < maxmax[post_area][post_type]){
			maxmax[post_area][post_type] = max_Isyn[post_nrn];	
		}
	}		
	
	for(int post_area = 0;post_area<N_areas;post_area++){
		for(int post_type = 0; post_type< Ntypes; post_type ++){
			
			if( glu_only && !glutamatergic[post_type] ) continue;
			
			//int br = Soma[post_nrn];
			
			for( int layer = 0; layer < Layers; layer ++ ){
				if( SYN[post_area][post_type][layer].size() > 0 ){  
					// synapses in that layer.
					
					if( maxmax[post_area][post_type] < limit ){
						// Scale down all incoming pathways.
						for( int pre_area = 0; pre_area < N_areas; pre_area++){
							for( int pre_type = 0; pre_type < MAX_TYPES; pre_type++){
								// Only adjust if we are in window for learning for each neural pathway.
								// (If higher order area is being scaled down before the lower areas have finished adjusting, the 
								//  higher order area may scale down too much and stop responding to the reduced activity in lower area
								//  once the lower area crystalizes after homeostasis.)  JLM
								if( simtime >= SYN[post_area][post_type][layer][pre_area][pre_type].lstarttimems && 
								   simtime < SYN[post_area][post_type][layer][pre_area][pre_type].lendtimems ){ 
									// Leave inhib neuron input alone.  We are trying to prevent excessive exc input.  If anything we should scale the inhib up.
									float target_sum = SYN[post_area][post_type][layer][pre_area][pre_type].init_max_sum;
									if( glutamatergic[pre_type] && target_sum > 0.0 ){
										target_sum = target_sum*(limit/maxmax[post_area][post_type]);
										
										SYN[post_area][post_type][layer][pre_area][pre_type].init_max_sum = target_sum;
										cout << "adjusting target sum to " << target_sum << ". Max current was "<<maxmax[post_area][post_type] << ". [post_area][post_type][layer][pre_area][pre_type]" << 
										post_area<<" "<<post_type<<" "<<layer<<" "<<pre_area<<" "<<pre_type << endl; // DEBUG 
									}
								}
							}
						}
					}
				}
			}
		}
	}
	
	for(int i = 0; i < max_Isyn.size();i++){
		max_Isyn[i]=0.0f;
	}
}

//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////

void normalize_synaptic_strengths(bool glu_only=false, bool initialization = false)
{ 
	
	// For each connection between types, normalize incoming synaptic strength.
	// That means making all synapses of a given type within a layer sum to a constant.
	
	for(int post_nrn=myid;post_nrn<N;post_nrn+=NP){
		
		int post_type = Type[post_nrn];
		
		if( glu_only && !glutamatergic[post_type] ) continue;
		
		int br = Soma[post_nrn];
		
		int post_area = Area[post_nrn];
		
		// look at all compartments in this neuron.
		while ( br<N_br && Br_N[br] == post_nrn ){
			
			int layer = Br_L[br];
			//vector<float> sum(MAX_TYPES,0.0f);
			vector<vector<float> > sum(N_areas, vector<float>(MAX_TYPES,0.0f));
			
			int first_br_in_layer = br;
			
			// look at all compartments in this layer.
			while ( br<N_br && Br_N[br] == post_nrn && Br_L[br] == layer ){
				
				for (int j=0;j<N_pre[br]; j++){
					int pre_type = Type[pre[br][j]];
					int pre_area = Area[pre[br][j]];
					
					//if( pre_area != post_area ) pre_type = Ntypes; // Fix for cortico-cortico connections.
					
					sum[ pre_area ][ pre_type ] += *s_pre[br][j];
					
				}
				//if(myid==0) cerr<< br << "*";
				
				br++;
			}
			
			// Now that we have the sums of all connection pathways in the layer, normalize them.
			
			br = first_br_in_layer;
			
			// predivide for efficiency.
			for( int pre_area = 0; pre_area < N_areas; pre_area++){
				for( int pre_type = 0; pre_type < MAX_TYPES; pre_type++){
					float target_sum = SYN[post_area][post_type][layer][pre_area][pre_type].init_max_sum;
					// only normalize total synaptic weight to the target value if simtime is larger than learning start time
					//  -- turn off synaptic weight until learning begins  - 01/13/2011  -- Y. Chen
					if (SYN[post_area][post_type][layer][pre_area][pre_type].lstarttimems>0) {
						if( simtime < SYN[post_area][post_type][layer][pre_area][pre_type].lstarttimems) {
							target_sum = 0.000001f;
						}
						sum[pre_area][pre_type] = target_sum/(sum[pre_area][pre_type] + 0.000001f);	
					}
					else {
						sum[pre_area][pre_type] = target_sum/(sum[pre_area][pre_type] + 0.000001f);
					}
				}
			}
			
			// look at all compartments in this layer
			while ( br<N_br && Br_N[br] == post_nrn && Br_L[br] == layer ){
				
				for (int j=0;j<N_pre[br]; j++){
					int pre_type = Type[pre[br][j]];
					int pre_area = Area[pre[br][j]];
					
					//if( pre_area != post_area ) pre_type = Ntypes; // Fix for cortico-cortico connections.
					
					*s_pre[br][j] *= sum[pre_area][pre_type];
					
					if ( initialization ){
						*s_pre[br][j] += *s_pre[br][j]*2*(random01()-.5)*SYN[post_area][post_type][layer][pre_area][pre_type].init_noise_pct;
					}
				}
				
				br++;
				
			}
			
			
			
		}
		
	}
	
	
	
	
}



//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////

void initialize_synaptic_strengths()
{ 
	int br, i, j, del;
	short post_type;
	
	
	for (br=0; br<N_br; br++)
	{
		
		post_type = Type[Br_N[br]];
		int post_neuron_number = Br_N[br];
		int post_area = Area[post_neuron_number];
		int layer = Br_L[br] ;
		float post_area_width_mm=area_data[ AreaNames[post_area] ].width_mm;
		
		
		for (j=0;j<N_pre[br]; j++)
		{
			int pre_neuron_number = pre[br][j];
			int pre_type = Type[pre_neuron_number];
			int pre_area = Area[pre_neuron_number];
			
			float d2,d_mm,sigma;
			
			connection conn= SYN[post_area][post_type][layer][pre_area][pre_type];
			sigma = conn.sigma_mm;		
			
			// if areas are different, project the neuron's axons onto the target area.
			
			float tx,ty;
			float tz=0.0f;
			float not_used;
			
			if( pre_area != post_area ){
				// Global connection.  Translate pre neuron's coord into local coordinates of target area.
				translate_pre_to_post_coord(pre_neuron_number, post_neuron_number, tx, ty, not_used);
			}
			else{
				// Areas are the same; no coordinate conversion necessary so use pre's coordinates as is; also no between-area delays.
				tx=x[pre_neuron_number][0];
				ty=x[pre_neuron_number][1];
			}
			
#ifdef WRAPAROUND
			d2 = (min(fabs(x[post_neuron_number][0]-tx),post_area_width_mm-fabs(x[post_neuron_number][0]-tx))*min(fabs(x[post_neuron_number][0]-tx),post_area_width_mm-fabs(x[post_neuron_number][0]-tx)))
			+(min(fabs(x[post_neuron_number][1]- ty),post_area_width_mm-fabs(x[post_neuron_number][1]-ty))*min(fabs(x[post_neuron_number][1]-ty),post_area_width_mm-fabs(x[post_neuron_number][1]-ty)))
			+(min(fabs(x[post_neuron_number][2]- tz),post_area_width_mm-fabs(x[post_neuron_number][2]-tz))*min(fabs(x[post_neuron_number][2]-tz),post_area_width_mm-fabs(x[post_neuron_number][2]-tz)));
#else
			d2 = (x[post_neuron_number][0]-tx)*(x[post_neuron_number][0]-tx) + (x[post_neuron_number][1]-ty)*(x[post_neuron_number][1]-ty) + (x[post_neuron_number][2]-tz)*(x[post_neuron_number][2]-tz);
#endif
			
			float dmin = conn.min_radius_mm;		
			float dmax = conn.max_radius_mm;		
			d_mm = sqrt(d2);
			if( dmin > .0001 ){
				// Surround type.  Calculate distance from ring.  Max of Gaussian connection strength is on the donut.
				float center_radius_mm = (dmin+dmax)*0.5f;
				d_mm = center_radius_mm - d_mm;
			}
			*s_pre[br][j]=exp(-(d_mm*d_mm)/(2*sigma*sigma));  // normal curve with max =1.0 at 0.0;
			*sd_pre[br][j]=0.0;
			*sds_pre[br][j]=0.0;
		}
		
	}	
	
	normalize_synaptic_strengths(false, true);
}


void reset_neuron( int nrn )
{
	int type = Type[nrn];
	int br = Soma[nrn];
	
	while(br>=0){
		
		Br_I_injected[br]=0.0;
		Br_AMPA[br]=0.0;
		Br_NMDA[br]=0.0;
		Br_NMDAVI[br]=0.0;
		Br_GABAa[br]=0.0;
		Br_GABAb[br]=0.0;
		Br_GABAc[br]=0.0;
		Br_V[br]=-60.0;
		//Br_U[br]=100.0*random01();  // NOTE: because there is rand # here, you cannot expect spike equivalence between reset_dynamics() and this method unless you call reset_neuron(nrn) in exctactly the order nrn=0:N 
		Br_U[br]=50.0;
		
		for (int j=0;j<sp_memory;j++) 
			dendritic_spikes[br][j]=-LT_plast;// fired a long time ago
		
		br=NextCompartment(br);
	 }

	for (int j=0;j<sp_memory;j++) 
			somatic_spikes[nrn][j]=-LT_plast;// fired a long time ago

	STSP_s[nrn][0]=1.0;
	STSP_s[nrn][1]=1.0;
		
	Background_AMPA[nrn]=Background_AMPA_for_type[type];
	Background_GABA[nrn]=Background_GABA_for_type[type];
	
	//note that unlike reset_dynamics(area,type) below, there's no clearing of the spike buffer here. so if you use this function in reset_interact directly instead of via any of the reset_dynamics(...) calls, you have to handle that manually
}


void reset_dynamics( group g )
{
	for (int nrn=g.first_neuron_id; nrn<=g.last_neuron_id; nrn++)
	{
		//if ( (nrn==g.first_neuron_id) && (myid == 0) )
		//	cerr << "Resetting group" << g.area_name+g.cell_name << " (" << g.first_neuron_id << "," << g.last_neuron_id <<")" << endl; 

		reset_neuron(nrn);
	}
		
	// it's ok if we do this multiple times on different reset calls, it's not much overhead and we don't do it often
	//DA=0.0; //RGM: don't want to reset		
	LR = 0.0; //RGM: Added
	B_firings=0;		// the first spike in the que
	N_firings=0;		// the last spike in the que
	firings_index[0]=-D;	// put a dummy spike at -D for simulation efficiency 
	firings_time[0]=0;	// timing of the dummy spike 
		
}

void reset_dynamics(int area, int type)
{
	if ( (area>=0) && (area<N_areas) && (type>=0) && (type<Ntypes) )
	{
		for(int nrn=itype[area][type]; nrn<itype[area][type+1]; nrn++)
		{
			reset_neuron(nrn);
		}
		
		// it's ok if we do this multiple times on different reset calls, it's not much overhead and we don't do it often
		//DA=0.0;	//RGM update	
		LR = 0.0; //RGM: Added
		B_firings=0;		// the first spike in the que
		N_firings=0;		// the last spike in the que
		firings_index[0]=-D;	// put a dummy spike at -D for simulation efficiency 
		firings_time[0]=0;	// timing of the dummy spike
		
		return;
	}
	
	cerr << "Ooops! Area " << area << "or Type " << type << "out of bounds in reset_dynamics(area,type)" << endl;
	exit(1);
}

void reset_dynamics(int area)
{
	for(int i=0; i<Ntypes; i++)
		reset_dynamics(area, i);	
}

void reset_dynamics()
{
	int br, i, j, del;
	short type, postclass;
	
	for (br=0; br<N_br; br++)
	{
		Br_I_injected[br]=0.0;
		Br_AMPA[br]=0.0;
		Br_NMDA[br]=0.0;
		Br_NMDAVI[br]=0.0;
		Br_GABAa[br]=0.0;
		Br_GABAb[br]=0.0;
		// Don't reset Br_GABAc
		Br_V[br]=-60.0;
		//Br_U[br]=100.0*random01(); 
		Br_U[br]=50.0;
	}
	
	for (i=0;i<N;i++)
		for (j=0;j<sp_memory;j++) 
			somatic_spikes[i][j]=-LT_plast;// fired a long time ago
	
	for (i=0;i<Br;i++)
		for (j=0;j<sp_memory;j++) 
			dendritic_spikes[i][j]=-LT_plast;// fired a long time ago
	
	for (j=0;j<STSP_types;j++) { // 		j=0 (no plasticity) should never be used in the code
		STSP[j][0] = 1;
		for (i=1;i<ST_plast; i++)
			STSP[j][i] = STSP[j][i-1]*(1.0 - 1.0/STSP_tau[j]);
	}
	
	
	for (i=0;i<N;i++) {
		STSP_s[i][0]=1.0;
		STSP_s[i][1]=1.0;
		int type = Type[i];
		
		Background_AMPA[i]=Background_AMPA_for_type[type];
		Background_GABA[i]=Background_GABA_for_type[type];
	}	
	
	
	//DA=0.0; //RGM update
	LR = 0.0; //RGM: Added
	
	B_firings=0;		// the first spike in the que
	N_firings=0;		// the last spike in the que
	firings_index[0]=-D;	// put a dummy spike at -D for simulation efficiency 
	firings_time[0]=0;	// timing of the dummy spike */
	
}


void initiate_dynamics()
{ 
	int br, i, j, del;
	short type, postclass;
	
	for (br=0; br<N_br; br++)
	{
		Br_I_injected[br]=0.0;
		Br_AMPA[br]=0.0;
		Br_NMDA[br]=0.0;
		Br_NMDAVI[br]=0.0;
		Br_GABAa[br]=0.0;
		Br_GABAb[br]=0.0;
		Br_GABAc[br]=0.0;
		Br_V[br]=-60.0;
		Br_U[br]=100.0*random01();
		
	}
	
	initialize_synaptic_strengths();
	
	for (i=0;i<N;i++)
		for (j=0;j<sp_memory;j++) 
			somatic_spikes[i][j]=-LT_plast;// fired a long time ago
	
	for (i=0;i<Br;i++)
		for (j=0;j<sp_memory;j++) 
			dendritic_spikes[i][j]=-LT_plast;// fired a long time ago
	
	
    LTD[0]=Aminus;
	LTD_sds[0] = Aminus_sds;
    for (i=1;i<LT_plast;i++)
	{
		LTD[i]=exp(-1.0/tau_minus)*LTD[i-1]; // making the LTD part exp(-dt/taum)
		LTD_sds[i] = exp(-1.0/tau_minus_sds)*LTD_sds[i-1];
	}
	
    LTP[0]=Aplus;
	LTP_sds[0] = Aplus_sds;
    for (i=1;i<LT_plast;i++)
	{
		LTP[i]=exp(-1.0/tau_plus)*LTP[i-1]; // making the LTP part exp(-dt/taup)
		LTP_sds[i] = exp(-1.0/tau_plus_sds)*LTP_sds[i-1];
	}
	
	for (j=0;j<STSP_types;j++) { // 		j=0 (no plasticity) should never be used in the code
		STSP[j][0] = 1;
		for (i=1;i<ST_plast; i++)
			STSP[j][i] = STSP[j][i-1]*(1.0 - 1.0/STSP_tau[j]);
	}
	
	
	for (i=0;i<N;i++) {
		STSP_s[i][0]=1.0;
		STSP_s[i][1]=1.0;
		int type = Type[i];
		
		Background_AMPA[i]=Background_AMPA_for_type[type];
		Background_GABA[i]=Background_GABA_for_type[type];
	}	
	
	DA=0.00; // RGM debug 0.0;
	LR = 0.0; //RGM: Added
	B_firings=0;		// the first spike in the que
	N_firings=0;		// the last spike in the que
	firings_index[0]=-D;	// put a dummy spike at -D for simulation efficiency 
	firings_time[0]=0;	// timing of the dummy spike
	
}


char fname_bra[30],fname_neu[30],fname_syn[30], fname_dyn[30];

void initialize()
{
	int	i,z;
	short type;
	
	// Dynamically allocate memory for s_pre, sd_pre, sds_pre, pre, & pre_delay
	// because their size Br * M is too large for the stack.
	s_pre = (float ***)malloc(Br*sizeof(float **));
	sd_pre = (float ***)malloc(Br*sizeof(float **));
	sds_pre = (float ***)malloc(Br*sizeof(float **));
	pre = (int **)malloc(Br*sizeof(int *));
	pre_delay = (unsigned char **)malloc(Br*sizeof(unsigned char *));
	
	if ( (s_pre == NULL) | (sd_pre == NULL) | (sds_pre == NULL) | (pre == NULL) | (pre_delay == NULL) )
	{
		fprintf(stderr,"Cannot allocate memory, processor myid=%d", myid);
		exit(0);
	}
	
	for (z = 0; z < Br; z++)
	{
		s_pre[z] = (float **)malloc(M*sizeof(float *));
		sd_pre[z] = (float **)malloc(M*sizeof(float *));
		sds_pre[z] = (float **)malloc(M*sizeof(float *));
		pre[z] = (int *)malloc(M*sizeof(int));
		pre_delay[z] = (unsigned char *)malloc(M*sizeof(unsigned char));
		
		if ( (s_pre[z] == NULL) | (sd_pre[z] == NULL) | (sds_pre[z] == NULL) | (pre[z] == NULL) | (pre_delay[z] == NULL) )
		{
			fprintf(stderr,"Cannot allocate memory, processor myid=%d", myid);
			exit(0);
		}
	}
	
	
	
	/*
	 itype[0]=0;	// the beginning of the first type
	 for (type=1; type<Ntypes; type++)
	 itype[type]=itype[type-1]+(int) floor(N*Types[type-1]/100.0);
	 itype[Ntypes]=N;
	 
	 for (type=0; type<Ntypes; type++) 
	 {	// for each neuronal type
	 
	 for (i=itype[type];i<itype[type+1];i++) 
	 { // for each neuron
	 Type[i]=type;
	 Soma[i]=-1;		// may not have a soma, if outside this subnetwork 
	 }
	 }
	 */
	
	// RAND_SEED is a hash define from the compile script (2nd argument)
	//srand(RAND_SEED); // reset the seed. All nodes need to generate exactly the same anatomy
	
	gsl_rng_set (rand_generator, RAND_SEED); // GSL
	
#ifdef new_somas
	allocate_somas_on_cortex();
	fprintf(stderr,"after allocate.");
	//sprintf(fname_neu,"%d-%d.neu",myid, N);
	//save_allocated_somas(fname_neu);
#else
	//load_allocated_somas(fname_neu);
#endif
	// N is now set.
	
	sprintf(fname_bra,"%d-%d.bra",myid, N);
	sprintf(fname_syn,"%d-%d.syn",myid, N);
	sprintf(fname_dyn,"%d-%d.dyn",myid, N);
	
	// RAND_SEED is a hash define from the compile script (2nd argument)
	//srand(RAND_SEED+myid+1);
	
	gsl_rng_set (rand_generator, RAND_SEED+myid+1); // GSL
	
#ifdef new_branches
	allocate_branches();
	save_allocated_branches(fname_bra);
#else
	load_allocated_branches(fname_bra);
#endif
	
	
	
	// RAND_SEED is a hash define from the compile script (2nd argument)
	//srand(RAND_SEED+myid);
	
	gsl_rng_set (rand_generator, RAND_SEED+myid); // GSL
	
#ifdef new_synapses
	if( myid == 0 ){ 
		for(int i = 0; i < N_areas; i++){
			write_connection_probabilities(AreaNames[i]); // For information only; not really necessary.
		}
	}
	allocate_synapses();
#else
	
#ifdef FindGlobalAxons
	allocate_synapses(); // need to establish S_local[br]
#endif
	
	load_allocated_synapses(fname_syn);
#endif
	
	
	// RAND_SEED is a hash define from the compile script (2nd argument)
	//srand(RAND_SEED+myid);
	
	gsl_rng_set (rand_generator, RAND_SEED+myid); // GSL
	
	initiate_dynamics();
	
#ifndef start_dynamics  // if not defined, then continue stimulation from checkpoint.
	load_dynamics(fname_dyn);
	
#endif
	
#ifdef new_synapses
	save_allocated_synapses(fname_syn);
#endif
	
	char fname[128];
	for (int area = 0; area < N_areas; area++)
	{
		for (type=0; type<Ntypes; type++) 
		{
			sprintf(fname,"voltages/%s%s",AreaNames[area],Names[type].c_str());
			remove(fname);  
		}
	}
	
#ifdef VARIABLE_TRANSMISSION
	neuron_count = N/NP + 1; // + 1 instead of ceiling, we don't want to truncate and miss xmitting something
	spike_buffer_out = (int *)malloc(sizeof(int)*neuron_count);
	spike_buffer_in = (int *)malloc(sizeof(int)*neuron_count*NP);
	
	if ( (spike_buffer_in == NULL) | (spike_buffer_out == NULL) )
	{
		cerr << "Memory allocation failed, processor myid=" << myid << endl;
		exit(0);
	}
	
	displacements[0] = 0;
	for (int i = 1; i < NP; i++)
		displacements[i] = displacements[i-1] + neuron_count;
#endif
#ifdef BITPACKED_TRANSMISSION
	nbytes_packed_out = N/NP/8 + 1; // 8 bits per byte, +1 also to avoid truncation
	nbytes_packed_in = nbytes_packed_out * NP;
	index_buffer = (bool *)malloc(sizeof(bool)*N);
	byte_buffer_out = (unsigned char *)malloc(sizeof(unsigned char)*nbytes_packed_out);
	byte_buffer_in = (unsigned char *)malloc(sizeof(unsigned char)*nbytes_packed_in);
	
	if ( (index_buffer == NULL) | (byte_buffer_out == NULL) | (byte_buffer_in == NULL) )
	{
		cerr << "Memory allocation failed, processor myid=" << myid << endl;
		exit(0);
	}
	
	for (int neu = 0; neu < N; neu++)
		index_buffer[neu] = false;
	
#endif
}



FILE *fvol[MAX_AREAS][MAX_TYPES];
FILE *debugvoltagefile;

void PrepareWritingVoltages ()
{	
	int type, nrn;
	char fname[128];
	// prepare to write voltage data
	for (int area = 0; area < N_areas; area++)
		for (type=0; type<Ntypes; type++) 
		{
#ifdef WRITEVOLTAGES1CELL			
			nrn = (itype[area][type]+itype[area][type+1])/2;	// one of the fist neuron of this type
			nrn = nrn - (nrn%NP);		// the first neuron on myid=0 of this type
			if( nrn < itype[area][type] ) nrn+=NP;
			
			sprintf(fname,"voltages/%s%s-%d",AreaNames[area],Names[type].c_str(),nrn);
#else
			nrn = itype[area][type];	// one of the fist neuron of this type
			nrn = nrn - (nrn%NP) + myid;		// the first neuron on myid=0 of this type
			if( nrn < itype[area][type] ) nrn+=NP;
			
			sprintf(fname,"voltages/%s%s-%d",AreaNames[area],Names[type].c_str(),myid);
			
#endif
			fvol[area][type] = fopen(fname,"w");
		}
	sprintf(fname,"voltages/debugvoltagefile-%d",myid);
	debugvoltagefile=fopen(fname,"w");
}

void WriteVoltages()
{
	
	int type, nrn, br;
	// save a complete copy of voltage distribtution for each neuronal type
	for (int area = 0; area < N_areas; area++)
		for (type=0; type<Ntypes; type++)
		{
			
#ifdef WRITEVOLTAGES1CELL			
			nrn = (itype[area][type]+itype[area][type+1])/2;	// one of the fist neuron of this type
			nrn = nrn - (nrn%NP);		// the first neuron on myid=0 of this type
			if( nrn < itype[area][type] ) nrn+=NP;
			br=Soma[nrn];
			
			
			while (br >=0 )
			{
				fprintf(fvol[area][type], "%d %d %3.3f %3.3f %3.3f %3.3f %3.3f %3.3f %3.3f   ",nrn, N_pre[br], Br_V[br], Br_U[br],Br_Isyn[br], Br_NMDA[br],Br_AMPA[br], Br_GABAa[br],Br_GABAb[br]);
				br = MoveAlongApicalDendrite(br);
				
			}
			fprintf(fvol[area][type], "\n");
#else
			// Write the soma values for all neurons on this processor.  This will slow things down.
			nrn = itype[area][type];	// one of the fist neuron of this type
			int procid_of_nrn=nrn%NP;
			nrn = nrn - procid_of_nrn + myid;		// the first neuron on myid=0 of this type
			if(nrn < itype[area][type]) nrn += NP;
			
			for( ; nrn < itype[area][type+1]; nrn += NP ){
				
				br=Soma[nrn];
				fprintf(fvol[area][type], "%d %d %3.3f %3.3f %3.3f %3.3f %3.3f %3.3f %3.3f   ",nrn, N_pre[br],Br_V[br], Br_U[br], Br_Isyn[br], Br_AMPA[br], Br_NMDA[br], Br_GABAa[br],Br_GABAb[br]);
				//fprintf(fvol[area][type], "%d %3.3f %3.3f ", nrn, Br_GABAa[br], Br_AMPA[br] );
				//fprintf(fvol[area][type], "%d %3.3f %3.3f ", nrn, Br_V[br], Br_U[br] );
				
			}
			fprintf(fvol[area][type], "\n");
#endif			
			
			
		}
	
	// DEBUG
	/*
	 nrn = 1050;
	 br=Soma[nrn];
	 while (br >=0 )
	 {
	 fprintf(debugvoltagefile, "%d %d %3.3f %3.3f %3.3f %3.3f %3.3f %3.3f %3.3f   ",nrn, N_pre[br], Br_V[br], Br_U[br],Br_Isyn[br], Br_NMDA[br],Br_AMPA[br], Br_GABAa[br],Br_GABAb[br]);
	 
	 br = NextCompartment(br); // returns -1 when no more compartments in the present neuron.
	 }
	 fprintf(debugvoltagefile, "\n");
	 */
}

void FinishWritingVoltages()
{
	int type, nrn;
	char fname[30], fname1[30];
	// close and rename files
	for (int area = 0; area < N_areas; area++)
		for (type=0; type<Ntypes; type++) 
		{
			fclose(fvol[area][type]); 
			
#ifdef WRITEVOLTAGES1CELL			
			nrn = (itype[area][type]+itype[area][type+1])/2;	// one of the fist neuron of this type
			nrn = nrn - (nrn%NP);		// the first neuron on myid=0 of this type
			if( nrn < itype[area][type] ) nrn+=NP;
			
			sprintf(fname,"voltages/%s%s-%d",AreaNames[area],Names[type].c_str(),nrn);
			sprintf(fname1,"voltages/%s%s-%d.dat",AreaNames[area],Names[type].c_str(),nrn);
#else
			nrn = itype[area][type];	// one of the fist neuron of this type
			nrn = nrn - (nrn%NP) + myid;		// the first neuron on myid=0 of this type
			if( nrn < itype[area][type] ) nrn+=NP;
			
			sprintf(fname,"voltages/%s%s-%d",AreaNames[area],Names[type].c_str(), myid);
			sprintf(fname1,"voltages/%s%s-%d.dat",AreaNames[area],Names[type].c_str(), myid);
#endif
			remove(fname1);
			rename(fname,fname1);  
		}
	fclose(debugvoltagefile); 
}


void WriteSynapticDistribution()
{
	int i,j, k, del, bin;
	int br;
	
#define histl 1000
	int hist[histl];
	
	float exex, exin;
	int	cex, cin;
	
	cex=0;
	cin=0;
	
	exex=0;
	exin=0;
	
	for (i=0;i<histl;i++) {
		hist[i]=0;
	}
	
	
	for (br=0;br<N_br;br++) 
	{
		for (j=0; j< N_pre[br]; j++) 
		{
			if (glutamatergic[ Type[ pre[br][j] ] ]==1) 
				if (glutamatergic[ Type[ Br_N[br] ]] ==1)
				{
					bin = floor(histl* (*s_pre[br][j])/50 );
					if (bin>=histl) bin=histl-1;
					hist[bin]++;
					exex += *s_pre[br][j];
					cex++;
				}
				else 
				{
					exin += *s_pre[br][j];
					cin++;
				}
		}
	}
	
	
	/*
	 for (i=0;i<N;i++)
	 if (glutamatergic[ Type[ i ] ] == 1) // yes, excitatory neuron
	 for (del=0;del<D;del++)
	 for (j=0;j<delays[i][del];j++) {
	 bin = floor(histl*s[i][del][j]/sm);
	 if (bin>=histl) bin=histl-1;
	 hist[bin]++;
	 }
	 */
	
	FILE *fs;
	fs = fopen("synhist.dat","a");
	for (i=0;i<histl;i++) 
		fprintf(fs, "%d ", hist[i]);
	fprintf(fs,"\n");
	
	fclose(fs);
	
	
	fs = fopen("synstr.dat","a");
	fprintf(fs, "%5.5f %5.5f\n", exex/cex, exin/cin);
	fclose(fs);
}


/*
 
 //*****************************************************************
 void homeostasis(){
 
 float scale_factor = 1.0f/(float)homeostasis_update_seconds;
 bool scale = false;
 float tolerance_ratio = 0.2; // 20% tolerance.
 
 //          					     nb1	p23	b23	 nb23 ss4(L4) ss4(L23) p4 b4	nb4	p5(L23)	p5(L56)	b5	nb5 p6(L4)	p6(L54)	b6	nb6	TCs	TCn	TIs	TIn	TRN
 
 vector<float> target_rate(Ntypes,1.0);
 
 
 if(myid==0) cout << "**************** APPLYING INTRINSIC AND SYNAPTIC HOMEOSTASIS **********************";
 
 // Homeostasis of background conductances to simulate intrinsic plasticity
 
 if(myid==0)cout << "(background_ampa, background_gaba)";
 
 for(int post_nrn=myid;post_nrn<N;post_nrn+=NP){
 if( (float)firing_rate[ post_nrn ]*scale_factor < target_rate[ Type[post_nrn] ] *(1.0 - tolerance_ratio) ){
 
 // Too low. scale up.
 Background_AMPA[post_nrn]=min(30.0, Background_AMPA[post_nrn]+0.25);
 Background_GABA[post_nrn]=min(30.0, Background_GABA[post_nrn]+0.25*1.5);
 }
 else if( (float)firing_rate[ post_nrn ]*scale_factor > target_rate[ Type[post_nrn] ] *(1.0 + tolerance_ratio ) ){
 
 // Too high.
 Background_AMPA[post_nrn]=max(0.0, Background_AMPA[post_nrn]-0.25);
 Background_GABA[post_nrn]=max(0.0, Background_GABA[post_nrn]-0.25*1.5);
 }
 if(myid==0)cout << '(' << Background_AMPA[post_nrn] << '.' << Background_GABA[post_nrn] << ") ";
 
 }
 if(myid==0)cout << endl;
 
 
 
 for(int br = 0; br<Br; br++){
 int post_nrn = Br_N[br];
 
 int posttype = Type[ post_nrn ];
 int postclass = Neuronal_Class[posttype];
 
 float exc_weight_scale = 1.0f, inh_weight_scale=1.0f;
 
 scale = false;
 if( (float)firing_rate[ post_nrn ]*scale_factor < target_rate[ Type[post_nrn] ] *(1.0 - tolerance_ratio) ){
 
 // Too low. scale up.
 exc_weight_scale = 1.05; // 5% increase every 10 seconds.
 inh_weight_scale = .95;
 scale  = true;
 }
 else if( (float)firing_rate[ post_nrn ]*scale_factor > target_rate[ Type[post_nrn] ] *(1.0 + tolerance_ratio ) ){
 
 // Too high.
 exc_weight_scale = .95;
 inh_weight_scale = 1.05;
 scale  = true;
 }
 
 // Scale it.
 if( scale ){
 //if (glutamatergic[Type[post_nrn]]==1){ // postsynaptic neuron is exc 
 for(int j = 0; j < N_pre[br]; j++){
 
 int pre_nrn = pre[br][j];
 int preclass = Neuronal_Class[ Type[pre_nrn] ];
 
 if (glutamatergic[Type[pre_nrn]]==1){ // presynaptic neuron is exc 
 *s_pre[br][j] *= exc_weight_scale;
 }
 else{
 *s_pre[br][j] *= inh_weight_scale;									
 }
 // clip
 if( *s_pre[br][j] > sm[preclass][postclass] ) *s_pre[br][j] = sm[preclass][postclass];
 if( *s_pre[br][j] < sm_min ) *s_pre[br][j] = sm_min;
 }
 //}
 }
 
 }
 
 // Clear the firing rate counters.
 for (int i=myid;i<N;i=i+NP) {
 if(myid==0)cout << setw(7) << fixed << setprecision(3) << firing_rate[i]*scale_factor << " ";
 firing_rate[i] = 0.0f;
 }
 if(myid==0)cout << endl;
 
 
 
 }
 
 */





//*************************************************************************************************
//*************************************************************************************************
// recreates the allocation of neurons on a processor, so we can know exactly how many neurons are on a given processor.
int Number_Neurons_On (int proc, int total_procs, int total_neurons)
{
	int n_temp, n_diff, n_per;
	n_per = ( (int) (total_neurons/total_procs) );
	n_temp = total_procs * n_per;
	n_diff = total_neurons - n_temp;
	
	if (proc < n_diff)
	{
		return (n_per + 1);
	}
	else
	{
		return n_per;
	}
}
//*************************************************************************************************
//*************************************************************************************************
void AddSpike(int nrn)
{
	int	pr, sp, nsp, i, k, tpassed;
	short	type, preclass, s_var, stsp_type;
	
	type = Type[nrn];
	
	frate[ type ]++;
	
	// For homeostasis: JLM 6-10-09
	firing_rate[nrn]+= 1.0;
	
	// Short-Term Synaptic Plasticity; 
	// While this spike is delivered to the postsynaptic targets, use the following 
	// value of the STSP_s variable to scale the synaptic strength
	
	preclass = Neuronal_Class[type];
	
	tpassed = simtime - somatic_spikes[nrn][0];  // current time - the last recorded spike, not counting this one
	if (tpassed >= ST_plast)
	{
		STSP_s[nrn][0] = 1.0;
		STSP_s[nrn][1] = 1.0;
	}
	else 
	{
		for (s_var=0; s_var<2; s_var++) 
		{
			stsp_type = STSP_s_type[s_var][preclass];
			if (stsp_type > 0) // i.e., there is some kind of plasticity
			{
				STSP_s[nrn][s_var] = 1 + (STSP_s[nrn][s_var]*STSP_p[stsp_type]-1)*STSP[stsp_type][tpassed];
				if (STSP_s[nrn][s_var] > 20)
					STSP_s[nrn][s_var] = 20; // manual cut-off
			}
		}
	}
	
	
	N_firings = (N_firings+1)%N_firings_max;
	if (N_firings == B_firings) 
	{
		if( myid == 0) fprintf(stderr,"UnpackSpikes(): Too many spikes at t=%d. Ignoring some.", t);
	}
	
	
	firings_time[N_firings] = t;		// the time of firing of the neuron
	firings_index[N_firings] = nrn;	// the identity of the fired neuron
	
	
	for (k=1;k<sp_memory;k++) somatic_spikes[nrn][k]=somatic_spikes[nrn][k-1];
	somatic_spikes[nrn][0] = simtime;
	
}
//*************************************************************************************************
//*************************************************************************************************
#ifdef FIXED_TRANSMISSION
void	UnpackSpikes(void)
{	
	int	pr, sp, nsp, nrn, frate;

#ifdef DISPLAY_WINDOW
  	int temp = 0;
  	N_spikes_for_display = 0;
#endif

#ifndef DEMO
	// write spikes
	FILE *fs, *fr;
	if(myid == 0 )
		fs = fopen("spikes.datt","a");
	if(myid == 0 )
		fr = fopen("frate.dat","a");
#endif
	
	frate=0;
	
	for( int i = 0; i < N; i++)
	{
		spiked[i]=false;
	}
	
	for (pr=0;pr<NP;pr++)
	{
		nsp = spikes_indeces_in[pr];	// the number of spikes from processor pr
#ifdef DISPLAY_WINDOW
		N_spikes_for_display = N_spikes_for_display + nsp;
#endif
		for (sp=0;sp<nsp;sp++) 
		{
			nrn=spikes_in[pr*max_transmit+sp];// index of fired neuron
			AddSpike(nrn);
			spiked[nrn]=true;
			frate++;
#ifndef DEMO
			if (myid == 0) 
				fprintf(fs, "%d  %d\n", simtime, nrn);
#endif
				
#ifdef DISPLAY_WINDOW
			sp_buffer[temp] = nrn;
			temp++;
#endif
		}
	}
	
#ifndef DEMO
	if(myid == 0) 
		fclose(fs);
	if(myid == 0 )
	{
		fprintf(fr, "%d\n", frate);
		fclose(fr);
	}
#endif
}
#endif
//*************************************************************************************************
//*************************************************************************************************
#ifdef VARIABLE_TRANSMISSION
void	VariableUnpackSpikes(void)
{	
	int	pr, sp, nsp, nrn, frate;

#ifdef DISPLAY_WINDOW
  	int temp = 0;
  	N_spikes_for_display = 0;
#endif	

#ifndef DEMO
	// write spikes
	FILE *fs, *fr;
	if(myid == 0 )
		fs = fopen("spikes.datt","a");
	if(myid == 0 )
		fr = fopen("frate.dat","a");
#endif	

	frate=0;
	
	for( int i = 0; i < N; i++)
	{
		spiked[i]=false;
	}
	
	for (pr=0;pr<NP;pr++)
	{
		nsp = spike_count_in[pr];	// the number of spikes from processor pr
#ifdef DISPLAY_WINDOW
		N_spikes_for_display = N_spikes_for_display + nsp;
#endif
		for (sp=0;sp<nsp;sp++) 
		{
			nrn = spike_buffer_in[ displacements[pr] + sp];	// index of fired neuron
			AddSpike(nrn);
			spiked[nrn] = true;
			frate++;
#ifndef DEMO
			if (myid == 0) 
				fprintf(fs, "%d  %d\n", simtime, nrn);
#endif
				
#ifdef DISPLAY_WINDOW
			sp_buffer[temp] = nrn;
		  	temp++;
#endif
		}
	}
	
#ifndef DEMO
	if(myid == 0) 
		fclose(fs);
	if(myid == 0 )
	{
		fprintf(fr, "%d\n", frate);
		fclose(fr);
	}
#endif
}
#endif
//*************************************************************************************************
//*************************************************************************************************
#ifdef BITPACKED_TRANSMISSION
void BitPackSpikes (void)
{
	int nrn;
	int pbyte = 0; 	// which byte we're on right now
	int pbit = 0; 		// number of bits we've processed
	unsigned char mask;
	
	memset(byte_buffer_out, 0x00, nbytes_packed_out);
	for (nrn = myid; nrn < N; nrn = nrn + NP) 
	{
		switch (pbit % 8)
		{
			case 0:
				mask=0x01;
				break;
			case 1:
				mask=0x02;
				break;
			case 2:
				mask=0x04;
				break;
			case 3:
				mask=0x08;
				break;
			case 4:
				mask=0x10;
				break;
			case 5:
				mask=0x20;
				break;
			case 6:
				mask=0x40;
				break;
			case 7:
				mask=0x80;
				break;
			default:
				fprintf(stderr,"Logic error in pack_spikes! pbit=%x\n", pbit);
				exit(-1);
		}
		
		if (index_buffer[nrn])
		{
			byte_buffer_out[pbyte] = byte_buffer_out[pbyte] | mask;
			
		}
		
		if (mask == 0x80) // reached end of this byte
			pbyte++;
		
		pbit++;	
	}
}
void BitUnpackSpikes (void)
{
	int nrn, frate;
	bool nrn_spiked;
	
	int pbyte = 0; 	// which byte we're on right now
	int pbit = 0; 		// number of bits we've processed from a proc
	unsigned char mask;
	
	FILE *fs, *fr;
	
#ifdef DISPLAY_WINDOW
	N_spikes_for_display = 0;
#endif

	frate = 0;
	
	for( int i = 0; i < N; i++)
	{
		spiked[i] = false;
	}
	
	if(myid == 0)
		fs = fopen("spikes.datt","a");
	if(myid == 0)
		fr = fopen("frate.dat","a");
		
	for (int proc = 0; proc < NP; proc++)
	{		
		nrn = proc;
		pbit = 0;
		pbyte = proc * nbytes_packed_out;
		
		while (pbit < Number_Neurons_On(proc, NP, N)) 
		{
			switch( pbit % 8)
			{
				case 0:
					mask=0x01;
					break;
				case 1:
					mask=0x02;
					break;
				case 2:
					mask=0x04;
					break;
				case 3:
					mask=0x08;
					break;
				case 4:
					mask=0x10;
					break;
				case 5:
					mask=0x20;
					break;
				case 6:
					mask=0x40;
					break;
				case 7:
					mask=0x80;
					break;
				default:
					fprintf(stderr,"Logic error in unpack_spikes! pbit=%x\n", pbit);
					exit(-1);
			}
			
			nrn_spiked = (bool) (byte_buffer_in[pbyte] & mask);
						
			if (nrn_spiked)
			{
				AddSpike(nrn);
				spiked[nrn] = true;
				frate++;
				
				if (myid == 0) 
					fprintf(fs, "%d  %d\n", simtime, nrn);
					
#ifdef DISPLAY_WINDOW
				sp_buffer[N_spikes_for_display] = nrn;
				N_spikes_for_display++;
#endif
			}

			if (proc == myid)
			{
			    	if (nrn_spiked && !index_buffer[nrn])
			    	{
					cerr << "Unpack Error: false false, on processor myid=" << myid << endl;
					exit(0);
				}
			    	if (!nrn_spiked && index_buffer[nrn])
			   	{
					cerr << "Unpack Error: false true, on processor myid=" << myid << endl;
					exit(0);
			   	}
			}
			
			if (mask == 0x80) // reached end of this byte
				pbyte++;
			
			index_buffer[nrn] = false;
			nrn = nrn + NP;
			pbit++;
			
		} // end of while we haven't processed all the spikes on this proc loop
		
	} // end of for all procs loop
	
	if(myid == 0) 
		fclose(fs);
	
	if(myid == 0) 
	{
		fprintf(fr, "%d\n", frate);
		fclose(fr);
	}
}
#endif
//*************************************************************************************************
//*************************************************************************************************




// Increase initial strength of TC cells (type 17[TCs] and 18 [TCn])
// Usage: 	Make_Syn_Stronger(17); Make_Syn_Stronger(18);
/*
 void Make_Syn_Stronger(int type)
 {
 int i, del, j;
 for( int a = 0;a<N_areas; a++)
 for (i=itype[a][type]; i<itype[a][type+1]; i++){
 int pre_class = Neuronal_Class[ Type[ i ] ];
 for (del=0; del<D;del++)
 for (j=0;j<delays[i][del];j++)
 s[i][del][j] = sm[pre_class][Neuronal_Class[   Type[Br_N[post[i][del][j]]]   ]]; 
 
 }
 }
 */


///////////////////////////  Lesion //////////////////////////
#ifdef Lesion
const	float Lesion_Center[3]={-0.3670,    0.6023,   -0.2568};
const	float	Lesion_Rad = 0.1;
short	Lesion_branch[Br]; // 0 - lesioned, 1 - OK

void Make_Lesion() // identify all neurons that belong to the lesion radius
{
	
	int br, i, k;
	float r;
	
	for (br=0; br<N_br; br++)
	{
		i=Br_N[br];
		r=0;
		for (k=0;k<3;k++)
			r += (x[i][k]-Lesion_Center[k])*(x[i][k]-Lesion_Center[k]);
		
		if (r<Lesion_Rad*Lesion_Rad)
			Lesion_branch[br]=0;
		else
			Lesion_branch[br]=1;
	}
} 

#endif




/*
 const int N_Electrodes = 25; // The number of electrodes
 
 const float Electrode_Center[][3] = { 
 -0.1284,    0.6591,    0.7545, //primary visual cortex (V1)
 };
 
 
 const float Electrode_rad = 0.08;
 
 
 //0="nb1","p23","b23","nb23","ss4(L4)","ss4(L23)","p4","b4","nb4",
 //9="p5(L23)","p5(L56)","b5","nb5","p6(L4)","p6(L56)","b6","nb6",
 //17="TCs","TCn","TIs","TIn",21="TRN";
 
 // Types of neurons that will be stimulated
 //const int Stim_Electrode_Ntypes = 3;	// ss4, ss4, TCs
 //const int Stim_Electrode_itypes[Stim_Electrode_Ntypes] = {4,5,17}; 
 //const int Stim_Electrode_Ntypes = 1;	// TCs
 //const int Stim_Electrode_itypes[Stim_Electrode_Ntypes] = {17};
 
 const int Stim_Electrode_Ntypes = 2;	// p4, TCs
 const int Stim_Electrode_itypes[Stim_Electrode_Ntypes] = {6, 17};
 
 // Types of neurons from which LFP's will be recorded
 const int Rec_Electrode_Ntypes = 12;	
 const int Rec_Electrode_itypes[Rec_Electrode_Ntypes] = {1, 2, 4, 5, 6, 7, 9, 10, 11, 13, 14, 15};
 
 
 const int Electrode_Max = (NMAX/N_Electrodes/NP);	// maximal number of neurons at each electrode 
 int	Stim_Electrode_N[N_Electrodes];						// the actual number of neurons at each electrode
 int 	Stim_Electrode[N_Electrodes][Electrode_Max];	// the indeces of neurons
 int	Rec_Electrode_N[N_Electrodes];						// the actual number of neurons at each electrode
 int 	Rec_Electrode[N_Electrodes][Electrode_Max];	// the indeces of neurons
 
 
 void Set_Electrodes() // identify all neurons that belong to the electrodes
 {
 
 int elec, type;
 int i, k;
 float r;
 
 
 for (elec=0;elec<N_Electrodes; elec++) 
 {
 Stim_Electrode_N[elec]=0;
 
 for (type = 0; type<Stim_Electrode_Ntypes; type++)
 for (i=itype[Stim_Electrode_itypes[type]]; i<itype[Stim_Electrode_itypes[type]+1]; i++)
 if (Soma[i]>=0) // i.e., within this subnetwork 	
 {
 r=0;
 for (k=0;k<3;k++)
 r += (x[i][k]-Electrode_Center[elec][k])*(x[i][k]-Electrode_Center[elec][k]);
 
 if (r<Electrode_rad*Electrode_rad)
 {
 if (Stim_Electrode_N[elec] < Electrode_Max)
 {
 Stim_Electrode[elec][Stim_Electrode_N[elec]++] = i;
 }
 }
 }	
 }
 
 for (elec=0;elec<N_Electrodes; elec++) 
 {
 Rec_Electrode_N[elec]=0;
 
 for (type = 0; type<Rec_Electrode_Ntypes; type++)
 for (i=itype[Rec_Electrode_itypes[type]]; i<itype[Rec_Electrode_itypes[type]+1]; i++)
 if (Soma[i]>=0) // i.e., within this subnetwork 	
 {
 r=0;
 for (k=0;k<3;k++)
 r += (x[i][k]-Electrode_Center[elec][k])*(x[i][k]-Electrode_Center[elec][k]);
 
 if (r<Electrode_rad*Electrode_rad)
 {
 
 if (Rec_Electrode_N[elec] < Electrode_Max)
 {
 Rec_Electrode[elec][Rec_Electrode_N[elec]++] = i;
 }
 }
 }	
 }
 remove("Electrodes.dat");
 
 
 } 
 
 float EEG_Electrodes[N_Electrodes];
 float EEG_Electrodes_in[NP*N_Electrodes];
 
 void Process_Electrodes()
 {
 int elec, i, br, pr;
 float V;
 
 for (elec=0;elec<N_Electrodes; elec++) 
 {
 EEG_Electrodes[elec] = 0.0;
 
 for (i=0; i<Rec_Electrode_N[elec]; i++)
 {
 
 br = Soma[ Rec_Electrode[elec][i] ];
 
 V = Br_V[br];
 while (br >=0 )
 {
 EEG_Electrodes[elec] += Br_V[br]-V;
 V=Br_V[br];
 br = MoveAlongApicalDendrite(br);
 }
 }
 
 if (Rec_Electrode_N[elec]>0)
 EEG_Electrodes[elec] = EEG_Electrodes[elec]/Rec_Electrode_N[elec];
 } 
 
 MPI_Allgather(EEG_Electrodes, N_Electrodes, MPI_FLOAT, EEG_Electrodes_in, N_Electrodes, MPI_INT,  MPI_COMM_WORLD); 
 
 if (myid != 0) return;
 
 FILE *fs;
 fs = fopen("Electrodes.dat", "a");
 
 for (elec=0;elec<N_Electrodes; elec++)
 {
 EEG_Electrodes[elec]=0.0;
 for (pr=0;pr<NP;pr++) 
 EEG_Electrodes[elec] += EEG_Electrodes_in[pr*N_Electrodes+elec];
 fprintf(fs,"%3.3f  ", EEG_Electrodes[elec]/NP);	
 } 	
 fprintf(fs,"\n");	
 fclose(fs);
 }
 
 */


void Set_Electrodes(){
	
	// Read the input cells from the inputs.dat file.
	
	FILE	*fs;
	fs = fopen("microstimulation_list.dat", "r");
	
	N_Inputs = 0;	
	if( fs!= NULL ){
		
		int neuron;
		while ( (N_Inputs<NMAX/2) && ( EOF != fscanf(fs,"%d ",&neuron)) )
		{
			if (Soma[neuron]>=0){ // i.e., within this subnetwork 	
				Inputs[N_Inputs++]=neuron;
				cout << neuron << " ";
			}
		}
	}
}





//////////////////   fMRI/BOLD signal ////////////


#ifdef fMRI

#define NNP 5000  			// the number of voxels for fMRI/BOLD signal
							// all voxels are some points (locations) in myid=0
							// any other position should be mapped to the target voxel
int	Voxel_target[N];		// the target voxel; i.e., this voxel data should be 
							// added to the target (nearlest) voxel in myid=0
							// only myid=0 voxels will be saved (still plenty)
int	Voxel_count[NNP];		// the count of how many voxels are averaged in this one


float BOLD_L[NNP]; 			// local BOLD signal
float BOLD_G[NNP];			// global BOLD signal
float A_GLUT[NNP];			// activity of glutamatergic synapses
float A_GABA[NNP];			// activity of gaba-ergic synapses

float	A_frate[Ntypes][NNP];// Firing rate of neurons of different types




void Setup_BOLD()
{
	float d, d_min;
	int	i, j, target;
	
	
	for (i=0;i<NNP; i++)
	{
		BOLD_L[i]=0.0;
		BOLD_G[i]=0.0;
		A_GLUT[i]=0.0;
		A_GABA[i]=0.0;
		
		for (j=0;j<Ntypes;j++)
			A_frate[j][i] = 0;
	}
	
	
	remove("bold/BOLD_L.dat");
	remove("bold/BOLD_G.dat");
	remove("bold/GLUT.dat");
	remove("bold/GABA.dat");
	
	char fname[30];
	for (j=0; j<Ntypes; j++)
	{
		sprintf(fname,"bold/vox_%s.dat",Names[j]);
		remove(fname);  
	}
	
	
	for (j=0;j<NNP;j++) 
		Voxel_count[j]=0;
	
	for (i=0;i<N;i++)
	{
		d_min = 10000;
		// find the nearest voxel in myid=0
		for (j=0; (j<N) & (j<NNP*NP) ; j=j+NP)
		{
			d=(x[i][0]-x[j][0])*(x[i][0]-x[j][0]) + (x[i][1]-x[j][1])*(x[i][1]-x[j][1]) + (x[i][2]-x[j][2])*(x[i][2]-x[j][2]);
			if (d<d_min)
			{
				d_min =  d;
				target = j;
			}
		}
		Voxel_target[i] = target/NP;
		Voxel_count[ Voxel_target[i] ]++;
	}
}



void Process_BOLD() // do this in the t-loop (every ms)
{
	for (int v=0; v<NNP; v++)
	{
		BOLD_L[ v ] *= (1-1.0/500);
		BOLD_G[ v ] *= (1-1.0/500);
		A_GLUT[ v ] *= (1-1.0/500);
		A_GABA[ v ] *= (1-1.0/500);
	}
}

float Voxels_in[NP*NNP];
float	Voxels_all[NNP];		// 


void Consolidate_voxels(float SIG[], char fname[30])
{
	
	int pr, bnrn, j;
	
	MPI_Allgather(SIG, NNP, MPI_FLOAT, Voxels_in, NNP, MPI_INT,  network_comm); 
	
	if (myid==0)
	{
		for (j=0;j<NNP; j++)
			Voxels_all[j]=0.0;
		
		for (j=0;j<NNP; j++)
		{
			for (pr=0;pr<NP;pr++)
			{
				// pr+j*NP is the global index of the voxel
				Voxels_all[j] += Voxels_in[pr*NNP+j]; 
			}
		}
		
		FILE *fs;
		fs = fopen(fname, "a");
		for (j=0;j<NNP; j++)
			fprintf(fs,"%3.3f  ", Voxels_all[j]/Voxel_count[j]);	
		
		fprintf(fs,"\n");	
		fclose(fs);
	}
	
}



void Save_BOLD()
{
	int type, i;
	
	Consolidate_voxels(BOLD_L, "bold/BOLD_L.dat");
	Consolidate_voxels(BOLD_G, "bold/BOLD_G.dat");
	Consolidate_voxels(A_GLUT, "bold/GLUT.dat");
	Consolidate_voxels(A_GABA, "bold/GABA.dat");
	
	
	char fname[30];
	for (type=0; type<Ntypes; type++)
	{
		sprintf(fname,"bold/vox_%s.dat",Names[type]);
		
		Consolidate_voxels(&(A_frate[type][0]), fname);
		for (i=0; i<NNP; i++)
			A_frate[type][i] = 0.0;
	}
	
}

#endif
////////////////////////////////////////////////

//**********************
// If there is an external group with the same area name and type as an internal
// group, then external spikes force the corresponding internal neuron to fire.
//
// WARNING: You should make sure that such local dummy groups do not have any synapses of their
// own which would couse them to fire by themselves.
class translator{
	
public:
	//*******************************
	// If the external neuron has no translation, returns -1
	// Otherwise converts from one to the other.
	int map_external_to_local_neuron_number(int external_neuron_id){
		
		if( external_to_internal_id_table.find(external_neuron_id) != external_to_internal_id_table.end() ){
			return external_to_internal_id_table[ external_neuron_id ];
		}
		else{
			return -1;
		}
		
	}
	
	translator(){
    
    MPI_Barrier(network_comm);
    
		internal_groups.load_groups("groups.dat");
		external_groups.load_groups("groups_external.dat");
		
		map<string,group> ext_groups = external_groups.get_map();
		
		max_spikes_per_msec=0;
		for( map<string,group>::iterator it = ext_groups.begin();it!=ext_groups.end();it++){
			if( internal_groups.exists( it->second.area_name+it->second.cell_name) ) {
				
				group ig = internal_groups[ it->second.area_name+it->second.cell_name ];
				
				if( ig.rows == it->second.rows && ig.cols == it->second.cols ){
					// All is well; put it in the lookup table.
					int internal_neuron_id= ig.first_neuron_id;
					for(int external_neuron_id = it->second.first_neuron_id ; external_neuron_id <= it->second.last_neuron_id; external_neuron_id++){ 
						external_to_internal_id_table[ external_neuron_id ] = internal_neuron_id;
						internal_neuron_id++;
						
						max_spikes_per_msec++; // keep track of externally transmitting neurons.
											   //cout << "<" << internal_neuron_id << ","<< external_neuron_id<<">"<<endl;
					}
				}
				else{
					
					cerr << "Warning: the external simulator has a matching group with a different size than the internal group with the same name"<< endl;
					cerr << "See groups_external.dat: " << ig.area_name << " " << ig.cell_name << "internal g:" << ig.rows << "," << ig.cols << " external g:" << it->second.rows << "," << it->second.cols << endl;
					cerr << "Either change the names if they should not match, or make the sizes match." << endl;					
					exit(1);					
				}
			}
		}
		max_spikes_per_msec = 10000;
		
	}
	
	// Assumes you can't transmit more spikes in 1 msec than the number of neurons mentioned in
	// groups_external.dat.
	int get_max_spikes_per_msec(){
		return max_spikes_per_msec;
	}
	
private:
	groups internal_groups;
	groups external_groups;
	int max_spikes_per_msec;
	map<int,int> external_to_internal_id_table;
};


float interp(float x1,float y1,float x2,float y2,float x){
	
	float y,m,b;
	b=x1;
	m=(y2-y1)/(x2-x1);
	return y1+m*(x-x1);
	
}


void run( int Tfin )
{
	
	int pr, i, j, k, del, arrived_ago, arrived_after;
	int soma, type, nrn, br, br_a, postbr, daughter, postnrn;
	short s_var, pretype, posttype, preclass, postclass, stsp_type;
	float stsp_s, PSc, syn_scale;
	float Idendr, Gdendr, Edendr, V, U, xx, NMDAgate, NMDAVIgate;
	float Gampa, Gnmda, Ggabaa, Ggabab, Eampa, Enmda, Egabaa, Egabab; 
	float Ggabac, Egabac;
	float E, F, g;
	float EEG, EEG_in[NP];
	FILE *fs;
	
	
	translator translatorx;
	int max_external_transmit=translatorx.get_max_spikes_per_msec();
	int* external_spikes = new int[max_external_transmit];
	
	
	int training_pts;
	double xmax=-1000000.0;
	double ymax=-1000000.0;
	
	float sm_min = 0.001;
	int pre_type, pre_area, post_type, post_area, layer;
	
	training_pts = read_training_array((char*)"a.dat",training_dim); // load training data from force_dv.dat;
	
	/*
	 for (int i=0;i<training_pts;i++) {
	 xmax=max(xmax,fabs(training_array[i][0]));
	 ymax=max(ymax,fabs(training_array[i][1]));
	 }
	 cout << xmax << " "  << ymax << " "  << training_pts << "\n";
	 
	 for (int i=0;i<training_pts;i++) {     // normalize training data to [0 1]
	 training_array[i][0]=training_array[i][0]/xmax/2.1+0.5;
	 training_array[i][1]=training_array[i][1]/ymax/2.1+0.5;
	 }
	 */
	
#ifdef LOG_GROUP
    FILE *filelog;
    //if(myid==0)
    char filelog_name[128];
    sprintf(filelog_name, "grouplog_%d.dat", myid);
    filelog=fopen(filelog_name,"w");

    int avg_count[LOG_COUNT];
    float avg_G[LOG_COUNT];
    float avg_E[LOG_COUNT];
    float avg_Gampa[LOG_COUNT]; 
    float avg_Gnmda[LOG_COUNT];
    float avg_Ggabaa[LOG_COUNT];
    float avg_Ggabab[LOG_COUNT];
    float avg_Ggabac[LOG_COUNT];
    float avg_Eampa[LOG_COUNT];
    float avg_Enmda[LOG_COUNT];
    float avg_Egabaa[LOG_COUNT];
    float avg_Egabab[LOG_COUNT];
    float avg_Egabac[LOG_COUNT];
    float avg_U[LOG_COUNT];
    float avg_V[LOG_COUNT];
    float avg_F[LOG_COUNT];
    float avg_Vr[LOG_COUNT];
    float avg_Vt[LOG_COUNT];
    
    float avg_syn[LOG_COUNT];
    int avg_nsyn[LOG_COUNT];
    
    int avg_temp[LOG_COUNT];
#endif
	
	fprintf(stderr,"Processor myid=%d finished\n",myid);
	if (myid == 0 )
	{
		remove("EEG.dat");
		remove("frates.dat");
		remove("synstr.dat");
		remove("synhist.dat");
		remove("frate.dat");
	}	
	
	// RAND_SEED is a hash define from the compile script (2nd argument)
	//srand(RAND_SEED);
	
	gsl_rng_set (rand_generator, RAND_SEED); // GSL
	
	Set_Electrodes();
	
#ifdef fMRI
	Setup_BOLD();
#endif
	
#ifdef Lesion
	Make_Lesion();
#endif
	
	
	
	
	
	/*  JLM: conflicts with homeostasis. 6-15-09
	 for (br=0; br<N_br; br++)
	 for (j=0;j<N_pre[br]; j++)
	 if (glutamatergic[Type[pre[br][j]]]==0) // presynaptic neuron is inh 
	 *s_pre[br][j]=sm[0];		// initial (and subsequent) inh values
	 */
	
	
	
	/*
	 
	 // Kill all the spikes, so that there is no ongoing activity
	 for (br=0; br<N_br; br++)
	 {
	 Br_AMPA[br]=0.0;
	 Br_NMDA[br]=0.0;
	 Br_GABAa[br]=0.0;
	 Br_GABAb[br]=0.0;
	 Br_V[br]=-60.0;
	 Br_U[br]=0;
	 }
	 B_firings=0;		// the first spike in the que
	 N_firings=0;		// the last spike in the que
	 firings_index[0]=-D;	// put a dummy spike at -D for simulation efficiency 
	 firings_time[0]=0;	// timing of the dummy spike
	 
	 */
	
	vector<float> max_Isyn(N,0.0);
	
	
	/*
	 for (i=0;i<N;i++)
	 if (glutamatergic[ Type[ i ] ] ==1) // yes, excitatory neuron
	 for (del=0;del<D;del++)
	 for (j=0;j<delays[i][del];j++)
	 {
	 posttype = Type[ Br_N[post[i][del][j]] ];
	 k = glutamatergic[ posttype ]; // 0-inh, 1-exc, exc -> ? synapse
	 if (k==0) // exc->inh
	 {
	 postclass = Neuronal_Class[ posttype ];
	 if (s[i][del][j]>sm[postclass]) 
	 s[i][del][j]=sm[postclass];
	 sd[i][del][j]=0;
	 }
	 }
	 */
	
	//if ( start_sec == 0) Make_Syn_Stronger(17);  // McKinstry; make sure the synapses are higher at the beginning of a simulation, not after restart.
	// main loop
	
	
#ifdef DISPLAY_WINDOW
	// signal display processor to continue & at the same time tell display processor
	// what the current time is.
	if (display_rank == 1)
	{
		MPI_Send (&start_sec, 1, MPI_INT, 0, 777, display_comm);
	}
#endif

#ifdef EXP3
			int tms;
			int trial_start_ms = g_motor_world::Instance()->trial_start_ms;
			int trial_duration = g_motor_world::Instance()->trial_duration;
			int lookat_A_ms = g_motor_world::Instance()->lookat_A;
			int lookat_B_ms = g_motor_world::Instance()->lookat_B;
			int lookaway_A_ms = g_motor_world::Instance()->lookaway_A;
			int lookaway_B_ms = g_motor_world::Instance()->lookaway_B;
#endif
	
	for (sec=start_sec;sec<Tfin;sec++)
	{
		
		//Make_Syn_Stronger(17); //Make_Syn_Stronger(18);
		
		// RAND_SEED is a hash define from the compile script (2nd argument)
		//srand(RAND_SEED+sec*NP+myid);
		
		gsl_rng_set (rand_generator, RAND_SEED+sec*NP+myid); // GSL
		
		if (myid == watchid) fprintf(stderr,"\nsec=%d ", sec);
#ifdef WRITEVOLTAGES
		if(sec % voltage_writing_frequency_seconds==0){
			
#ifdef WRITEVOLTAGES1CELL			
			if (myid == 0) PrepareWritingVoltages();
#else
			PrepareWritingVoltages();		
#endif
		}
#endif		
		
#ifndef DEMO
		if (myid == 0 )remove("spikes.datt");
#endif
		if (myid == 0) WriteSynapticDistribution();
		
		for (i=0;i<N;i++) 
			spiked[i]=false; // no spikes at the beginning.
		
		for (t=0;t<1000;t++)
		{

#ifdef APEX3
	if(myid == 0)
	{
		APEX_button_press = 0;
		changemode(1);
		if( kbhit() )
		{
			int keypress = getchar();
			int exit = false;
			std::cout << "Registered Keypress: " << keypress << endl;
			
			if(keypress == 32) 
			{
				std::cout << "PAUSE BUTTON PRESSED!!!!!!!!!!!!!!!!!" << std::endl;
				//pause simulation
				while( !exit ) 
				{ 
					if( kbhit() == 1 ) exit = true;
					std::cout << " . " << std::endl; sleep(1); 
				} 
			}
			else if( keypress == 59 )
			{
				std::cout << "APEX REWARD BUTTON PRESSED!!!!!!!!!!!!!!!!!!" << std::endl;
				APEX_button_press = 1;
			}
		}
		changemode(0);
	}	
	MPI_Bcast(&APEX_button_press, 1, MPI_INT, 0, network_comm);
	APEX_REWARD = APEX_button_press;
	
	if( APEX_REWARD )
	{
		std::cout << "REWARD SET!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
	}

#endif
      //if( myid == watchid) fprintf(stderr,".");
			simtime = sec*1000+t;
			
			N_spikes_out=0;
			
			//Process_Electrodes();
			EEG = 0.0;
			
#ifdef fMRI
			Process_BOLD();
#endif
			
			// Update the state of neurons and branches
			// Two ways this could in princple be done:
			// loop for neurons 1-N, for each neuron a loop through compartments
			// loop for each compartment 1-N_br. Let us do the latter.
			// ********************************************
			
#ifdef Stimulate
			
			//for (i=myid;i<N;i=i+NP) 
			//	exc_input[i] = 0.0f;
			for (i=0;i<N;i=i++) 
				exc_input[i] = 0.0f;
			
			int training_index= simtime/250;  //training data sampling rate is 100Hz, 10ms per point -- Y. Chen
			int cycle_repeat = training_index/training_pts;
			int data_index = training_index - cycle_repeat*training_pts;
			// if (i==0) 
			//cout << training_pts << " " << training_index << " " << cycle_repeat << " " << data_index << " "  << training_array[data_index][0] << " " << training_array[data_index][1] << "\n"; // debug for index
			
			//motor_world.interact(simtime, spiked, exc_input,training_array[data_index],training_array[data_index]);		
#ifdef EXP3
			g_motor_world::Instance()->interact(simtime, spiked, exc_input, training_array[data_index],training_array[data_index]);
#else
			motor_world.interact(simtime, spiked, exc_input, training_array[data_index],training_array[data_index]);
#endif
			
			//***********************************************
			// Read spikes from external simulator(s)
			int n_external_spikes;
			
			if(myid == 0)
			{
				//cout << "Just before MPI_Bcast" << endl;
				MPI_Bcast(&n_external_spikes, 1, MPI_INT, 0, exchange_comm);
        //RGM: DEBUG
        /*if( n_external_spikes == 0 ){
          cerr << "No External Spikes! =========================================" << endl << endl;
        }*/
				if( n_external_spikes>max_external_transmit ){
					cerr << "Too many spikes to send from external simulator!.  Ignoring extra spikes." << endl;
					n_external_spikes=max_external_transmit;
				}
				MPI_Bcast(external_spikes, n_external_spikes, MPI_INT, 0, exchange_comm);
				MPI_Barrier(exchange_comm);
				//MPI_Bcast(&n_external_spikes, 1, MPI_INT, 1, exchange_comm); // for synchronization purposes; don't know why MPI_Barrier didn't work.
				//MPI_Barrier(exchange_comm);  // Without this, the external simulator may get ahead creating buffer overflow.
				//fputs("received broadcast\n", stderr);
			}
			
#ifdef EXP3
			tms = sec*1000 + t;
			// no visual input during DELAY-A
			if( ((tms-trial_start_ms)%trial_duration >= lookaway_A_ms) && ((tms-trial_start_ms)%trial_duration < lookat_B_ms) )
			{
				n_external_spikes = 0;
			}
			// no visual input during DELAY-B
			if( (tms-trial_start_ms)%trial_duration >= lookaway_B_ms )
			{
				n_external_spikes = 0;
			}
#endif
			
			// Now convert the neuron numbers to numberes in my simulation and broadcast external spikes to all 
			MPI_Bcast(&n_external_spikes, 1, MPI_INT, 0, network_comm);
			MPI_Bcast(external_spikes, n_external_spikes, MPI_INT, 0, network_comm);
      MPI_Barrier(network_comm);//RGM
			// Pretend like all of these spikes were generated in my simulation.
			/*
			 cerr << "myid: " << myid;
			 for (i = 0; i < n_external_spikes; ++i)
			 fprintf(stderr, " (%d) %d,", i, external_spikes[i]);
			 */
			// Where do these neurons map to?
			for(i=0;i<n_external_spikes;i++){
				int neuron_id = translatorx.map_external_to_local_neuron_number(external_spikes[i]);
				
				//cerr << "myid: " << myid;
				//fprintf(stderr, " (%d) %d,", i, neuron_id);
				
				if( neuron_id >= 0 && Soma[ neuron_id] >=0 ){ // if this external spike is of interest to us && neuron is onboard, make dummy spike.
															  //cerr << "myid: " << myid;
															  //fprintf(stderr, " (%d) %d,", i, neuron_id);
					
					Br_V[ Soma[ neuron_id ] ]=100;  // Enough to force a spike.
				}
			}
			
			
			//***********************************************
			
			
			
			for (i=myid;i<N;i=i+NP){
				//	if( i>= 17640 && i<=17670 && (t%250 == 200) ){
				//		cout << i << "@" << myid << "-" << exc_input[i] << " ";	
				//	}
				//if( i>= 17640 && i<=17670 && (t%250 == 200) ) cout << endl;
				//if( exc_input[i] != 0.0f ){
				Br_I_injected[Soma[i]] = exc_input[i];
				//Br_NMDA[Soma[i]] = Br_NMDA[Soma[i]] + exc_input[i];
				//}
			}
			
			// reset_interact must be defined in world.h;  it uses the helper functions reset_XXXX that are defined above.  To maintain the original functionality, simply define
			// void reset_interact(int t) {
			//		if (t%pattern_duration_msec == 0)
			//			reset_dynamics();
			// }
			// otherwise you may control resetting in detail by using reset_dynamics(area) or reset_dynamics(area,type) or reset_neuron(nrn).
#ifdef EXP3
			g_motor_world::Instance()->reset_interact(simtime);
#else
			motor_world.reset_interact(simtime);
#endif

			
			if (((sec-start_sec)%60)<30) {  			// 30sec of stim, 30 sec of silence 
				
				
				// Stimulation
				
				if (sec >=1)
					//			if (t%100 == getrandom(100))
				{
					if (myid==0 && N_Inputs>0) fprintf(stderr," STIM ");
					
					for (i=0; i< N_Inputs; i++)
					{
						
						br = Soma[ Inputs[i] ];
						
						float stim_current;
						if (t%500<250){
							stim_current = 200.0f; // /(1-1.0/5.0); // compensate for decay.
						}
						else{
							stim_current = 0.0f;
						}
						
						Br_I_injected[br] = stim_current;
						
						Br_AMPA[br] = 0.0;
						Br_NMDA[br] = 0.0;
						Br_NMDAVI[br] = 0.0;
						
						Br_GABAa[br] = 0.0;
						Br_GABAb[br] = 0.0;
						Br_GABAc[br] = 0.0;
						
						
						//Br_AMPA[br] = Br_AMPA[br] + stim_current;
						//Br_NMDA[br] = Br_NMDA[br] + stim_current;
						
						//Br_V[br]=100;
					}
				}
				
				
			}
			
			
#endif
			
#ifdef LOG_GROUP
      for (int log_c=0; log_c < LOG_COUNT; log_c++) {
        avg_count[log_c] = 0;
        avg_G[log_c] = 0.0;
        avg_E[log_c] = 0.0;
        avg_Gampa[log_c] = 0.0; 
        avg_Gnmda[log_c] = 0.0;
        avg_Ggabaa[log_c] = 0.0;
        avg_Ggabab[log_c] = 0.0;
        avg_Ggabac[log_c] = 0.0;
        avg_Eampa[log_c] = 0.0;
        avg_Enmda[log_c] = 0.0;
        avg_Egabaa[log_c] = 0.0;
        avg_Egabab[log_c] = 0.0;
        avg_Egabac[log_c] = 0.0;
        avg_U[log_c] = 0.0;
        avg_V[log_c] = 0.0;
        avg_F[log_c] = 0.0;
        avg_Vr[log_c] = 0.0;
        avg_Vt[log_c] = 0.0;
        
        avg_syn[log_c] = 0.0;
        avg_nsyn[log_c] = 0;
        
        avg_temp[log_c] = 0;
      }
#endif
			
#define AchMod 0

			for (br=0;br<N_br;br++)
			{
				nrn = Br_N[br];
				type = Type[nrn];
				soma = Soma[nrn];
				syn_scale =1.0; // maybe decreased for dendrites of inhibitory neurons (below)
				
#ifdef LOG_GROUP
        V = Br_V[br];
				U = Br_U[br];
        
        Gampa=syn_scale*(Br_AMPA[br]+Background_AMPA[nrn]);
				Gnmda=syn_scale*(NMDAgate*Br_NMDA[br]+NMDAVIgate*Br_NMDAVI[br]);
				Ggabaa=syn_scale*(Br_GABAa[br]+Background_GABA[nrn]);
				Ggabab=syn_scale*(AchMod+Br_GABAb[br]);
				Ggabac=syn_scale*(Br_GABAc[br]);
				
				Eampa=syn_scale*E_glu*(Br_AMPA[br]+Background_AMPA[nrn]);
				Enmda=syn_scale*E_glu*(NMDAgate*Br_NMDA[br]+NMDAVIgate*Br_NMDAVI[br]);
				Egabaa=syn_scale*E_GABAa*(Br_GABAa[br]+Background_GABA[nrn]);
				Egabab=syn_scale*E_GABAb*(AchMod + Br_GABAb[br]);
				Egabac=syn_scale*E_GABAc*(Br_GABAc[br]);
        
        E =  - U + Edendr + (Eampa+Enmda+Egabaa+Egabab+Egabac) + Br_I_injected[br];
				g =  Gdendr + (Gampa+Gnmda+Ggabaa+Ggabab+Ggabac);        
				F = K[type]*(V-vr[type])*(V-vt[type]) ;
        
        post_area = Area[ Br_N[br] ];
        post_type = Type[ Br_N[br] ];
        for (int log_c=0; log_c<LOG_COUNT; log_c++) {
          if ( (simtime >= LOG_START_MS[log_c]) && (simtime < LOG_STOP_MS[log_c]) && (post_area == GROUP_AREAS[log_c]) && (post_type == GROUP_TYPES[log_c]) ) {
            avg_count[log_c]++;
            avg_G[log_c] += g;
            avg_E[log_c] += E;
            avg_Gampa[log_c] += Gampa; 
            avg_Gnmda[log_c] += Gnmda;
            avg_Ggabaa[log_c] += Ggabaa;
            avg_Ggabab[log_c] += Ggabab;
            avg_Ggabac[log_c] += Ggabac;
            avg_Eampa[log_c] += Eampa;
            avg_Enmda[log_c] += Enmda;
            avg_Egabaa[log_c] += Egabaa;
            avg_Egabab[log_c] += Egabab;
            avg_Egabac[log_c] += Egabac;
            avg_U[log_c] += U;
            avg_V[log_c] += V;
            avg_F[log_c] += F;
            avg_Vr[log_c] += vr[type];
            avg_Vt[log_c] += vt[type];
            
            for (j=0; j<N_pre[br]; j++) {
              if ( (Area[ pre[br][j] ] == 1) && (Type[ pre[br][j] ] == 18) ) { // pre: V1TCs
                avg_syn[log_c] += *s_pre[br][j];
                avg_nsyn[log_c]++;
              }
            }
          } else {
            avg_temp[log_c]++;
          }
        }
#endif

				if (br == soma ) 
				{
					
					
					//compute global EEG
					/*
					 if (glutamatergic[ Type[ nrn ]]==1) 
					 {
					 V = Br_V[soma];
					 br_a = MoveAlongApicalDendrite(soma);
					 while (br_a >=0 )
					 {
					 EEG += Br_V[br_a]-V;
					 V=Br_V[br_a];
					 br_a = MoveAlongApicalDendrite(br_a);
					 }
					 }
					 */
					
					if (simtime >= gabac_stop_msec[type]){
						Br_GABAc[soma] = 0.0; 
					}
					
					if (Br_V[soma] >= spike_peak_soma[type])
					{
						
						Br_V[soma] = spike_reset_c_soma[type]; // i.e., parameter c in the simple model
						Br_U[soma] += spike_reset_d_soma[type]; // i.e., parameter d in the simple model
						
						if ( (simtime >= gabac_start_msec[type]) && (simtime < gabac_stop_msec[type]) ){
							Br_GABAc[soma] += gabac_const[type]; // nS
						} 
						
						// post branch fired after all pre neurons
						for (j=0;j<N_pre[soma];j++) 
						{
							for (k=0;k<sp_memory;k++)
							{
                                arrived_ago = simtime-pre_delay[soma][j]-somatic_spikes[ pre[soma][j] ][k];
                                if ((arrived_ago >=0) && (arrived_ago < LT_plast) ) // find the most recent presynaptic spike, considering the delay
                                {
                                    *sd_pre[soma][j] += LTP[arrived_ago];
									
									// ST-STDP addition.
									pre_type	= Type[ pre[soma][j] ];
									pre_area	= Area[ pre[soma][j] ];
									post_type	= Type[ Br_N[soma] ];
									post_area	= Area[ Br_N[soma] ];
									layer		= Br_L[ soma ];
									
									switch (SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_type) {
										case NONE:
											break;
										case REGULAR:
											*sds_pre[soma][j] += LTP_sds[arrived_ago];
											break;
										case LINEAR_REDUCTION:
											if ( *sds_pre[soma][j] > 0)
												*sds_pre[soma][j] += (1 - (*sds_pre[soma][j]/SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max)) * LTP_sds[arrived_ago];
											else 
												*sds_pre[soma][j] += LTP_sds[arrived_ago];
											
											if ( *sds_pre[soma][j] > SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max )
												*sds_pre[soma][j] = SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max;
											break;
										case POLY_REDUCTION:
											if ( *sds_pre[soma][j] > 0)
												*sds_pre[soma][j] += (1 - 
																	  (*sds_pre[soma][j]/SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max)*
																	  (*sds_pre[soma][j]/SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max) )
																	  * LTP_sds[arrived_ago];
											else 
												*sds_pre[soma][j] += LTP_sds[arrived_ago];
											
											if ( *sds_pre[soma][j] > SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max )
												*sds_pre[soma][j] = SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max;
											break;
										default:
											fprintf(stderr,"Unknown ST-STDP type: %d",SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_type);
											exit(0);
											break;
									}
									
									break;
                                }
							}
						}
						
						
						
						// This is needed for LTD. We treat somatic compartment as a dendritic one for the sake of LTD
						for (k=1;k<sp_memory;k++) dendritic_spikes[soma][k]=dendritic_spikes[soma][k-1];
						dendritic_spikes[soma][0]=simtime;
						
#ifdef FIXED_TRANSMISSION
						if (N_spikes_out<max_transmit)
						{
							spikes_out[N_spikes_out++] = nrn; // to be transmitted to others (and self)
						}
						else
						{
							if( myid == watchid) 
								fprintf(stderr,"Too many spikes at t=%d. Ignoring the rest. ", t);
						}
#endif						
#ifdef VARIABLE_TRANSMISSION
						spike_buffer_out[N_spikes_out++] = nrn;
#endif
#ifdef BITPACKED_TRANSMISSION
						index_buffer[nrn] = true;
#endif
						
						// we skip the rest for br - it just fired a spike
						continue; 
					}
				}
				else
				{	// check for spike in dendritic compartments (soma was checked above)
					//if (glutamatergic[ type ] ==0)
					//	syn_scale =0.5;
					
					
					if (Br_V[br]>spike_peak_dend[type])
					{
						Br_V[br] = spike_reset_c_dend[type];
						Br_U[br] += spike_reset_d_dend[type]; 
						
						for (k=1;k<sp_memory;k++) dendritic_spikes[br][k]=dendritic_spikes[br][k-1];
						dendritic_spikes[br][0]=simtime;
						
						
						// post branch fired after all pre neurons
						for (j=0;j<N_pre[br];j++) {
							for (k=0;k<sp_memory;k++)
							{
                                arrived_ago = simtime-pre_delay[br][j]-somatic_spikes[ pre[br][j] ][k];
                                if ( (arrived_ago>=0) && (arrived_ago < LT_plast) ) // find the most recent presynaptic spike, considering the delay
                                {
									*sd_pre[br][j] += LTP[arrived_ago];
									
									// ST-STDP addition.
									pre_type	= Type[ pre[br][j] ];
									pre_area	= Area[ pre[br][j] ];
									post_type	= Type[ Br_N[br] ];
									post_area	= Area[ Br_N[br] ];
									layer		= Br_L[ br ];
									
									switch (SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_type) {
										case NONE:
											break;
										case REGULAR:
											*sds_pre[br][j] += LTP_sds[arrived_ago];
											break;
										case LINEAR_REDUCTION:
											if ( *sds_pre[br][j] > 0)
												*sds_pre[br][j] += (1.0 - (*sds_pre[br][j]/SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max)) * LTP_sds[arrived_ago];
											else 
												*sds_pre[br][j] += LTP_sds[arrived_ago];
											
											if ( *sds_pre[br][j] > SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max )
												*sds_pre[br][j] = SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max;
											break;
										case POLY_REDUCTION:
											if ( *sds_pre[br][j] > 0)
												*sds_pre[br][j] += (1.0 - 
																	(*sds_pre[br][j]/SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max)*
																	(*sds_pre[br][j]/SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max)
																	);
											else 
												*sds_pre[br][j] += LTP_sds[arrived_ago];
											
											if ( *sds_pre[br][j] > SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max )
												*sds_pre[br][j] = SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_max;
											break;
										default:
											fprintf(stderr,"Unknown ST-STDP type: %d",SYN[post_area][post_type][layer][pre_area][pre_type].ststdp_type);
											exit(0);
											break;
									}
									
									break;
                                }
							}
						}
						
						// we skip the rest for br - it just fired a spike
						continue; 
					}
				}
				
				
				
				//				if( myid == watchid) fprintf(stderr,"time %d, br=%d\n", t, br );
				
				
				V = Br_V[br];
				U = Br_U[br];
				xx = (-80-V)/60;
				xx = xx*xx;
				if (glutamatergic[type]==1)
					NMDAgate = xx/(1+xx);
				else
					NMDAgate = 1.0*xx/(1+xx);

				// This will simulate roughly Voltage Independent NMDA (VI) channels found between layer IV spiny stellate cells.
				xx = (-100-V)/60;
				xx = xx*xx;
				if (glutamatergic[type]==1)
					NMDAVIgate = xx/(1+xx);
				else
					NMDAVIgate = 1.0*xx/(1+xx);
				
				
				Br_AMPA[br] *= decay_ampa[type];
				Br_NMDA[br] *= decay_nmda[type];
				Br_NMDAVI[br] *= decay_nmda[type];
				Br_GABAa[br] *= decay_gabaa[type];
				Br_GABAb[br] *= decay_gabab[type];
				Br_GABAc[br] *= decay_gabac[type];
				
				// Dendritic currents coming from daughter compartments
				Gdendr = 0.0;
				Edendr = 0.0;
				
				daughter=Br_Up[br];
				float temp;
				temp = Gdendrites_Up[type];
				while (daughter >= 0)
				{
					//Idendr += Gdendrites_Up[type]*(V-Br_V[daughter]);
					//Idendr = Gdendr*V-Edendr;
					//if( Br_Sister[daughter] < 0 && Br_Up[daughter] >= 0) temp = Gdendrites_Up[type] *10.0;  // Make the apical trunk thicker by factor of 2 (conductance = pi*compartment_radius^2/(intracellular_resistance*compartment_length_microns).
					Gdendr += temp;
					Edendr += temp*Br_V[daughter];
					
					daughter = Br_Sister[daughter]; // go to the next daughter brahch.
				}
				
				// Dendritic current coming from mother compartment
				if (Br_Down[br] >=0) // there is a mother branch (i.e., br is not a soma)
				{
					//Idendr += Gdendrites_Down[type]*(V-Br_V[Br_Down[br]]);
					//Idendr = Gdendr*V-Edendr;
					Gdendr += Gdendrites_Down[type];
					Edendr += Gdendrites_Down[type]*Br_V[Br_Down[br]];
					
				}
				
				/* //Ach modulation turned off
				 #define AchModulation 0
				 float		AchMod;
				 if (glutamatergic[ Type[ nrn ] ] == 1)
				 AchMod = AchModulation;
				 else
				 AchMod = 0;
				 
				 AchMod *=  1 - 1.0/(1+3*(sec+0.001*t - start_sec) ); // smooth transition 
				 
				 
				 if (glutamatergic[ Type[ nrn ] ] == 1)
				 syn_scale *= 1 - 0.7*(1 - 1.0/(1+1*(sec+0.001*t - start_sec) ));
				 
				 */ 
				
				
				Idendr = Gdendr*V - Edendr;
				
				//Isyn=(V-E_glu)*(Br_AMPA[br]+NMDAgate*Br_NMDA[br])+(V-E_GABAa)*Br_GABAa[br]+(V-E_GABAb)*Br_GABAb[br]; // synaptic current
				//Isyn=V*(Gampa+Gnmda+Ggabaa+Ggabab) - (Eampa+Enmda+Egabaa+Egabab);
				Gampa=syn_scale*(Br_AMPA[br]+Background_AMPA[nrn]);   // Br_AMPA is conductance. 
				Gnmda=syn_scale*(NMDAgate*Br_NMDA[br]+NMDAVIgate*Br_NMDAVI[br]);
				Ggabaa=syn_scale*(Br_GABAa[br]+Background_GABA[nrn]);
				Ggabab=syn_scale*(AchMod+Br_GABAb[br]);
				Ggabac=syn_scale*(Br_GABAc[br]);
				
				
				Eampa=syn_scale*E_glu*(Br_AMPA[br]+Background_AMPA[nrn]);
				Enmda=syn_scale*E_glu*(NMDAgate*Br_NMDA[br]+NMDAVIgate*Br_NMDAVI[br]);
				Egabaa=syn_scale*E_GABAa*(Br_GABAa[br]+Background_GABA[nrn]);
				Egabab=syn_scale*E_GABAb*(AchMod + Br_GABAb[br]);
				Egabac=syn_scale*E_GABAc*(Br_GABAc[br]);
				Br_Isyn[br]=V*(Gampa+Gnmda+Ggabaa+Ggabab+Ggabac) - (Eampa+Enmda+Egabaa+Egabab+Egabac);
				
				
				
				//"Implicit" numerical scheme
				//V' 		= ( 3*(V+60)*(V+50) - U - Idendr - Isyn )/C;
				//V(t+1) = V(t) + tau* ( 3*(V(t)+60)*(V(t)+50) - U(t) - Gdendr*V(t+1)+Edendr - V(t+1)*(Gsyn)+(Esyn) )/C;
				//E =  - U(t) + Edendr + (Esyn) ;
				//F =  3*(V(t)+60)*(V(t)+50) ;
				//g =  Gdendr + (Gsyn) ;
				//V(t+1) = V(t) + tau*( F + E - V(t+1)*g  )/C
				//V(t+1)*(1+tau*g/C) = V(t)+tau*(F+E)/C;
				//V(t+1) = ( V(t) +tau*(F+E)/C )/(1+tau*g/C); 
				
				E =  - U + Edendr + (Eampa+Enmda+Egabaa+Egabab+Egabac) + Br_I_injected[br];
				g =  Gdendr + (Gampa+Gnmda+Ggabaa+Ggabab+Ggabac);
        
				//V' = ( 0.5*(V+60)*(V+50) - U - Idendr - Isyn )/20.0;
				F = K[type]*(V-vr[type])*(V-vt[type]) ;
				V = ( V + (F+E)/C[type] )/(1+g/C[type]);
				V=min(max(-90.0f,V),51.0f);  // Make up for integration errors.
				float old_V=Br_V[br]; // save for special cases below.
				Br_V[br] = V;
				Br_U[br] += a[type]*(b[type]*(V-vr[type])-U);
				if (Br_U[br] > U_upper_bound[type])
					Br_U[br] = U_upper_bound[type];
				
				
				// Wierd special cases that have to be inserted in the code.
				switch (type)
				{
					case 0:  // nb1 type (non-basket interneurons in Layer 1)
							 // LS: Late-spiking neurons
							 //V' = ( 0.3*(V+66)*(V+40) - U - Idendr - Isyn )/20.0;
						if (br != soma) // this is the soma compartment
						{
							Br_V[br] = old_V + ( -Idendr - Br_Isyn[br] )/250.0;
						}
						break;
						
					case 17: // TCs
							 //V' = ( 1.6*(V+60)*(V+50) - U - Idendr - Isyn )/200.0;
						if (V<-65)
							Br_U[br] = U + 0.01*(15*(V+65)-U);
						else
							Br_U[br] = U + 0.01*(-U);
						
						break;
					case 18: // TCn
							 //V' = ( 1.6*(V+60)*(V+50) - U - Idendr - Isyn )/200.0;
						if (V<-65)
							Br_U[br] = U + 0.01*(15*(V+65)-U);
						else
							Br_U[br] = U + 0.01*(-U);
						
						break;
						
						
					case 21: // TRN
							 //V' = ( 0.25*(V+65)*(V+45) - U - Idendr - Isyn )/40.0;
						if (V<-65)
							Br_U[br] = U + 0.015*(10*(V+65)-U);
						else
							Br_U[br] = U + 0.015*(2*(V+65)-U);
						break;
						//default:			
				}
				
			}
      
#ifdef LOG_GROUP
      for (int log_c=0; log_c<LOG_COUNT; log_c++) {
        //if ( (myid==0) && (avg_count[log_c]>0) && (simtime >= LOG_START_MS[log_c]) && (simtime < LOG_STOP_MS[log_c]) ) {
        if ( (avg_count[log_c]>0) && (simtime >= LOG_START_MS[log_c]) && (simtime < LOG_STOP_MS[log_c]) ) {
          avg_G[log_c] /= avg_count[log_c];
          avg_E[log_c] /= avg_count[log_c];
          avg_Gampa[log_c] /= avg_count[log_c]; 
          avg_Gnmda[log_c] /= avg_count[log_c];
          avg_Ggabaa[log_c] /= avg_count[log_c];
          avg_Ggabab[log_c] /= avg_count[log_c];
          avg_Ggabac[log_c] /= avg_count[log_c];
          avg_Eampa[log_c] /= avg_count[log_c];
          avg_Enmda[log_c] /= avg_count[log_c];
          avg_Egabaa[log_c] /= avg_count[log_c];
          avg_Egabab[log_c] /= avg_count[log_c];
          avg_Egabac[log_c] /= avg_count[log_c];
          avg_U[log_c] /= avg_count[log_c];
          avg_V[log_c] /= avg_count[log_c];
          avg_F[log_c] /= avg_count[log_c];
          avg_Vr[log_c] /= avg_count[log_c];
          avg_Vt[log_c] /= avg_count[log_c];
          if( avg_nsyn[log_c] > 0 ) {
            avg_syn[log_c] /= avg_nsyn[log_c];
          }

          fprintf(filelog,"t=%d, group=%s%s, count=%d, else=%d, N_br=%d, SynW=%f, V=%f, U=%f, F=%f, G=%f, E=%f, Vr=%f, Vt=%f, Gampa=%f, Gnmda=%f, Ggabaa=%f, Ggabab=%f, Ggabac=%f, Eampa=%f, Enmda=%f, Egabaa=%f, Egabab=%f, Egabac=%f\n", simtime, AreaNames[ GROUP_AREAS[log_c] ], Names[ GROUP_TYPES[log_c] ].c_str(), avg_count[log_c], avg_temp[log_c], N_br, avg_syn[log_c], avg_V[log_c], avg_U[log_c], avg_F[log_c], avg_G[log_c], avg_E[log_c], avg_Vr[log_c], avg_Vt[log_c], avg_Gampa[log_c], avg_Gnmda[log_c], avg_Ggabaa[log_c], avg_Ggabab[log_c], avg_Ggabac[log_c], avg_Eampa[log_c], avg_Enmda[log_c], avg_Egabaa[log_c], avg_Egabab[log_c], avg_Egabac[log_c]);
        }
      }
#endif
      
			// ********************************************
			// end of update of states of neurons
#ifdef WRITEVOLTAGES
			if(sec % voltage_writing_frequency_seconds==0){
#ifdef WRITEVOLTAGES1CELL
				if (myid == 0) WriteVoltages(); 
#else
				WriteVoltages(); 
#endif
			}
#endif
			//if( myid == 0) fprintf(stderr,"sec %d, t=%d, p=%s, wrote voltages\n", sec, t, processor_name );fflush(stderr);
			
			
			/*
			 // Exchange global EEG
			 MPI_Allgather(&EEG, 1, MPI_FLOAT, EEG_in, 1, MPI_FLOAT,  network_comm); 
			 
			 if( myid == watchid) {
			 
			 EEG=0;
			 for (pr=0;pr<NP; pr++)
			 EEG += EEG_in[pr];
			 
			 fs = fopen("EEG.dat", "a");
			 fprintf(fs, " %4.8f ", EEG/N);
			 fclose(fs);
			 }	
			 */
			
			
			
			
			
			// Now, trigger synaptic transmission for arrived spikes
			for (i=0;i<MINIs*M*N_br/1000;i++) // add minis
			{
				br = getrandom(N_br);
				
				if(N_pre[br]<1)continue;  // sometimes we have branches with no synapses near the fringes of an area.
				
				posttype = Type[ Br_N[br] ];
				postclass = Neuronal_Class[posttype];
				
				if (glutamatergic[ Type[ Br_N[br] ] ] == 1) 
				{
					j = getrandom(N_pre[br]);
					
					nrn = pre[br][j];
					pretype = Type[ nrn ];
					
					preclass = Neuronal_Class[pretype];
					stsp_type = STSP_Type[preclass][postclass];	// The type of synapse (from STP point of view)
					if (stsp_type == 0) // no short-term synaptic plasticity
						stsp_s = 1.0;
					else
					{
						s_var = STSP_s_var[stsp_type];					// whether we use the first or the second variable
						stsp_s = STSP_s[nrn][s_var];						// the scaling factor
					}
					
					// ?!?!?!
					// ST-STDP addition.
					pre_area	= Area[ nrn ];
					post_area	= Area[ Br_N[br] ];
					layer		= Br_L[ br ];
					
					PSc = stsp_s*(*s_pre[br][j])*(1.0 + SYN[post_area][posttype][layer][pre_area][pretype].ststdp_modulation*(*sds_pre[br][j]));	// sds addition for ST-STDP
					
					if ( PSc > 20*SYN[post_area][posttype][layer][pre_area][pretype].max )
						PSc = 20*SYN[post_area][posttype][layer][pre_area][pretype].max;
					if ( PSc < 0 )
						PSc = 0;
					
					if (glutamatergic[ pretype ]==1) {
						Br_AMPA[br] += PSc; // add rand to simulate minis
						Br_NMDA[br] += PSc; 
						// Mini's don't add anything to NMDAVI channels?
					} else {
						Br_GABAa[br] += PSc; // add rand to simulate minis
						Br_GABAb[br] += PSc; // add rand to simulate minis
					}
				}
			}
			
			
			

#ifdef FIXED_TRANSMISSION
			MPI_Allgather(spikes_out, max_transmit, MPI_INT, spikes_in, max_transmit, MPI_INT, network_comm); 
			MPI_Allgather(&N_spikes_out, 1, MPI_INT, spikes_indeces_in, 1, MPI_INT, network_comm); 
			
			if( N_spikes_out != spikes_indeces_in[myid] )
				fprintf(stderr,"myid=%d, sec %d, t=%d, Warning! AllGather returned #spikes different than I sent!\n", myid, sec, t );fflush(stderr);
			
			UnpackSpikes();
#endif
#ifdef VARIABLE_TRANSMISSION
			MPI_Allgather(&N_spikes_out, 1, MPI_INT, spike_count_in, 1, MPI_INT, network_comm); 
			MPI_Allgatherv(spike_buffer_out, N_spikes_out, MPI_INT, spike_buffer_in, spike_count_in, displacements, MPI_INT, network_comm); 
			
			if( N_spikes_out != spike_count_in[myid] )
				fprintf(stderr,"myid=%d, sec %d, t=%d, Warning! AllGather returned #spikes different than I sent!\n", myid, sec, t );fflush(stderr);
			
			VariableUnpackSpikes();
#endif
#ifdef BITPACKED_TRANSMISSION
			BitPackSpikes ();
			MPI_Allgather (byte_buffer_out, nbytes_packed_out, MPI_BYTE, byte_buffer_in, nbytes_packed_out, MPI_BYTE, network_comm);
			BitUnpackSpikes ();
#endif

			
			
#ifdef DISPLAY_WINDOW
			if (myid == 0)
			{
				MPI_Send (&N_spikes_for_display, 1, MPI_INT, 0, 700, display_comm);

				if (N_spikes_for_display > display_max_transmit)
			    		N_spikes_for_display = display_max_transmit;

			  	if (N_spikes_for_display > 0)
			    		MPI_Send (sp_buffer, N_spikes_for_display, MPI_INT, 0, 701, display_comm);
				
#ifdef EXP3
				g_motor_world::Instance()->apex_client.transmit_to_display();
#endif
			}
#endif
			
			// TESTING slow decrease in receptive field size as in SOM by slowly bringing inhibition online.
			float inhib_gain=1.0f;
			//			if( sec <= 25 ){
			//				inhib_gain = 0.0f; // interp(0.0f, 0.0f, (float)32, 1.0f,(float)sec);
			//			}
			//			else{
			//				inhib_gain = 1.0f;  // 2nd phase, no extra kick of DA.
			//			}			
			//			if(myid==0)cout << "???????????????????????? inhib_gain: " << inhib_gain << endl;
			
			
			// Deliver EPSCs and IPSCs here
			
			k=N_firings;		// start with the most recent spike and go backwards
			while ( (t-firings_time[k]<D) && (k != B_firings))  // i.e., consider only spikes sent at most D ms ago
			{
				
				del = t-firings_time[k];	// the delay from the spike
				int prenrn = firings_index[k];		// the pre-synaptic neuron
				pretype = Type[prenrn];			// the type of pre-neuron
				preclass = Neuronal_Class[pretype];
				
#ifdef fMRI
				A_frate[pretype][	Voxel_target[ prenrn ] ] += 1.0;
#endif
				
				for (j=0; j<delays[prenrn][del]; j++)
				{
					postbr=post[prenrn][del][j];	// postsynaptic branch
					postnrn = Br_N[postbr];		// postsynaptic neuron
					posttype = Type[postnrn];	// type of postsynaptic neuron
					postclass = Neuronal_Class[posttype];
					
					
					stsp_type = STSP_Type[preclass][postclass];	// The type of synapse (from STP point of view)
					if (stsp_type == 0) // no short-term synaptic plasticity
						stsp_s = 1.0;
					else
					{
						s_var = STSP_s_var[stsp_type];					// whether we use the first or the second variable
						stsp_s = STSP_s[prenrn][s_var];						// the scaling factor
					}
					
					// ?!?!?!
					// ST-STDP addition.
					pre_area	= Area[ prenrn ];
					post_area	= Area[ postnrn ];
					layer		= Br_L[ postbr ];
					
					PSc = stsp_s*s[prenrn][del][j]*(1 + SYN[post_area][posttype][layer][pre_area][pretype].ststdp_modulation*sds[prenrn][del][j]);	// sds addition for ST-STDP
					
					if ( PSc > 20*SYN[post_area][posttype][layer][pre_area][pretype].max )
						PSc = 20*SYN[post_area][posttype][layer][pre_area][pretype].max;
					if ( PSc < 0 )
						PSc = 0;
					
					if (glutamatergic[pretype]==1)
					{
						
						
#ifdef Lesion
						if (Lesion_branch[postbr] == 0) // i.e., it is lesioned of global connections
						{
							// Let's find out whether the synapse is local or global
							for (short jj=N_pre[postbr];jj>=S_local[postbr]; jj--)
								if (pre[postbr][jj]==prenrn)
									PSc=0.0; // no input
						}
#endif
						
						
#ifdef fMRI
						// Adding to the BOLD Signal (excitatory neurons)
						
						A_GLUT[ Voxel_target[Br_N[postbr]]  ] += PSc;
						
						short Bglobal=0;
						for (short jj=N_pre[postbr];jj>=S_local[postbr]; jj--)
							if (pre[postbr][jj]==prenrn)
							{ 
								Bglobal=1;
								BOLD_G[ Voxel_target[Br_N[postbr]] ] += PSc;
							}
						
						if ( Bglobal==0 )
							BOLD_L[ Voxel_target[Br_N[postbr]]  ] += PSc;
#endif						
						post_area= Area[postnrn];
						pre_area = Area[prenrn];
						layer = Br_L[postbr] ;
						
						float ampa_gain = SYN[post_area][posttype][layer][pre_area][pretype].ampa_gain;
						float nmda_gain = SYN[post_area][posttype][layer][pre_area][pretype].nmda_gain;
						float nmdaVi_gain = SYN[post_area][posttype][layer][pre_area][pretype].nmdaVi_gain;  
						   						
						Br_AMPA[postbr] += ampa_gain*PSc;
						
						// DEBUG
						//if( Br_N[postbr] == 1050 ){
						//		cout << "t=" << t << " pre:"<<prenrn << ":" << stsp_s << " " << s[prenrn][del][j] << " " << stsp_s*s[prenrn][del][j] << " " << Br_AMPA[postbr] << endl;
						//}
						
						Br_NMDA[postbr] += nmda_gain*PSc;  // Chen 2-12-10
						Br_NMDAVI[postbr] += nmdaVi_gain*PSc;  // Chen 2-12-10
						
						if (glutamatergic[ Type[ Br_N[postbr] ] ]==1) // if the post neurons is excitatory
						{
							// this spike is before postsynaptic spikes (always)
							arrived_after = simtime-dendritic_spikes[postbr][0];
							if (arrived_after<LT_plast)
							{
								sd[prenrn][del][j] -= LTD[ arrived_after ]; 
								
								// ST-STDP addition.
								if (SYN[post_area][posttype][layer][pre_area][pretype].ststdp_type != NONE) {
									sds[prenrn][del][j] -= LTD_sds[ arrived_after ];
								}
								
							}
							
							// DEBUG
							
							//if(prenrn >= 3871 && postnrn <=1155 && postbr==Soma[postnrn]){
							//if(prenrn == 3585 && postnrn ==2275 ){
							//	cout << "************ sd: " << sd[prenrn][del][j] << " arrived_after " << arrived_after << " " << flush;
							//}
							
							
						}
						
						
					}
					else
					{
						
#ifdef fMRI			// Inhibitory interneurons; always local BOLD
						BOLD_L[ Voxel_target[Br_N[postbr]]  ] += PSc;
						A_GABA[ Voxel_target[Br_N[postbr]]  ] += PSc;
						
#endif
						
						//PSc *= inhib_gain;  // trying to simulate a wider receptive field early, narrowing over time.
						
						// For inhib neurons these parameters are treated as exclusion window in which this pathway is ineffective.
						// This allows us to create epilepsy for a brief time which helps with map formation (we think.).
						
						post_area= Area[postnrn];
						
						pre_area = Area[prenrn];
						
						layer = Br_L[postbr] ;
						
						//if( simtime > SYN[post_area][posttype][layer][pre_area][pretype].lendtimems ||
						   
						//   simtime < SYN[post_area][posttype][layer][pre_area][pretype].lstarttimems ){ 
							
							float gabab_gain = SYN[post_area][posttype][layer][pre_area][pretype].gabab_gain;
							
							Br_GABAa[postbr] += PSc;
							Br_GABAb[postbr] += gabab_gain* PSc;  // McKinstry, 7-1-09.  According to Shepherd book, GABAb input conductance is weaker, but because
														   // reversal potential is -90, it has the same IPSP magnitude as GABAa currents.
						//}
						
						
					}
				}
				
				if (--k<0)
					k=N_firings_max-1; // The mod function, %,  does not go from <0 to >0 
				
			}
			B_firings=k;
			
			
			
			// Update excitatory synaptic weights every 10 ms
			//DA *= 0.995; // decay of dopamine, added the decay only in the equation below
			//cout << "Dopamine Decay Test, Dopamine term = " << DA << endl << endl;
			
			
#ifdef Allow_Plasticity
			
			
			// Find the maximum depolarizing current (I) input to each cell.
			if(simtime%1000>500) {
			 for (i=myid;i<N;i+=NP){
			
				int br_id = Soma[i];
				if( Br_Isyn[br_id] < max_Isyn[i] ){
					max_Isyn[i] = Br_Isyn[br_id];
				}
			 }
			}

			if( simtime % 1000  == 999  && sec >= homeostatic_scaling_adjustment_start_sec){
				bool glu_only=true;
				homeostatic_scaling_adjustment(max_Isyn,glu_only);
			}
			
			
			// ST-STDP (sds)
			if ( t%10 == 0 )
			{
				for (i = 0; i < N; i++)
				{
					if (glutamatergic[Type[i]] == 1)
					{
						for (del = 0; del < D; del++)
						{
							for (j = 0; j < delays[i][del]; j++)
							{
								sds[i][del][j] *= 0.998;
							}
						}
					}
				}
			}
			
			if (t%50==0)
			{
				if ( sec >= 0 ){
					for (i=0;i<N;i++){
						pre_type = Type[i];
						pre_area = Area[i];
						if( glutamatergic[ Type[ i ] ] == 1 ) // yes, excitatory neuron
							for( del=0;del<D;del++ )
								for( j=0;j<delays[i][del];j++ )
								{
									post_type = Type[ Br_N[post[i][del][j]] ];
									post_area = Area[ Br_N[post[i][del][j]] ];
									layer = Br_L[post[i][del][j]] ;
																		
									k = glutamatergic[ post_type ]; // 0-inh, 1-exc, exc -> ? synapse
									if( k==1 ) // exc->exc
									{
										/* Dopamine Code Addition (RGM) */
										if( simtime >= SYN[post_area][post_type][layer][pre_area][pre_type].DAstarttimems && 
										   simtime < SYN[post_area][post_type][layer][pre_area][pre_type].DAendtimems )
										{									
											//if(myid==0)cout << "???????????????????????? DA: " << DA << endl;
																						
											s[i][del][j]+=sd[i][del][j]*DA;
											
											if( s[i][del][j] > SYN[post_area][post_type][layer][pre_area][pre_type].max ) 
											{
												s[i][del][j]= SYN[post_area][post_type][layer][pre_area][pre_type].max;
											}
											if( s[i][del][j] < sm_min )
											{ 
												s[i][del][j]= sm_min;	
											}	
										}
										//else if normal plasticity
										else if( simtime >= SYN[post_area][post_type][layer][pre_area][pre_type].lstarttimems && 
										   simtime < SYN[post_area][post_type][layer][pre_area][pre_type].lendtimems )
										{ 										
											float t0=(float)SYN[post_area][post_type][layer][pre_area][pre_type].lstarttimems;
											float t1=(float)SYN[post_area][post_type][layer][pre_area][pre_type].lendtimems;
											
											
											//s[i][del][j]+=sd[i][del][j]*(0.25+DA);
											float lrate0=SYN[post_area][post_type][layer][pre_area][pre_type].lrate0;
											float lrateN=SYN[post_area][post_type][layer][pre_area][pre_type].lrateN;
											
											// TESTING slow reduction in learning rate as in SOM
											if(simtime<t0+(t1-t0)/2.0)
											{
												LR = interp(t0, lrate0, t0+(t1-t0)/2.0, lrateN,(float)simtime); 
											}
											else
											{
												LR = lrateN; 
											}									
											
											s[i][del][j]+=sd[i][del][j]*LR; //RGM: *DA
											// DEBUG
											//if( i==7400 && Br_N[post[i][del][j]]>3968 && Br_N[post[i][del][j]]<4050 ) cout << DA << " *************" << endl;
											//if(simtime>=5990 )cout << DA << " *************" << endl;
											
											if (s[i][del][j]> SYN[post_area][post_type][layer][pre_area][pre_type].max) 
												s[i][del][j]= SYN[post_area][post_type][layer][pre_area][pre_type].max;
											if (s[i][del][j]< sm_min) 
												s[i][del][j]= sm_min;
										} 
										//this term is used for depression or facilitation for the STDP rule depedending on the ordering of pre/post cells                                
										sd[i][del][j] *= sd_decay_rate;//0.951 // try decaying this value faster . . . could be trailing behind from the last stimulus
									}//end if k==1
								}
					}
					bool glut_only = false;
					//if (simtime%25000==0) 
					normalize_synaptic_strengths(glut_only);
					DA *= DA_decay_rate;
				} 
				else{  // give network a moment to settle down.
					
					for (i=0;i<N;i++)
						if (glutamatergic[ Type[ i ] ] ==1) // yes, excitatory neuron
							for (del=0;del<D;del++)
								for (j=0;j<delays[i][del];j++)
								{
									int posttype = Type[ Br_N[post[i][del][j]] ];
									k = glutamatergic[ posttype ]; // 0-inh, 1-exc, exc -> ? synapse
																   //					if (k==1) // exc->exc
									{
										sd[i][del][j]=0.0;
									}
								}
					
				}
				//				fprintf(stderr,"myid=%d finished plasticity\n",myid);
				
			}
			
#endif
			
#ifdef fMRI
			if (t%1000==0)
				Save_BOLD();
#endif
			
			
			
			
		} // end of t loop (1000ms = 1sec)
		
		
		
#ifdef WRITEVOLTAGES
		if(sec % voltage_writing_frequency_seconds==0){
			
#ifdef WRITEVOLTAGES1CELL
			if (myid == 0) FinishWritingVoltages();
#else
			FinishWritingVoltages();
#endif
		}
#endif
		
		// find the firing rate
		if( myid == watchid) fprintf(stderr,"\nFiring rate:\n ");
		
		if( myid == watchid) {
			fs = fopen("frates.dat", "a");
			for (type=0;type<Ntypes;type++) {
				fprintf(stderr,"[%s]=%3.3f, ", Names[type].c_str(), frate[type]/(N*Types[type]/100));
				fprintf(fs,"%3.3f ", frate[type]/(N*Types[type]/100) );
				frate[type]=0.0;
			}		
			fprintf(fs, "\n");
			fclose(fs);
		}	
		
		if( sec % homeostasis_update_seconds == homeostasis_update_seconds -1 ){
			//homeostasis();
			// DEBUG
			//sprintf(fname_syn,"%d-%dfinal.syn",myid, N);
			//save_allocated_synapses(fname_syn);
		}
		if( sec % save_synapses_period_seconds == save_synapses_period_seconds -1 ){
			sprintf(fname_syn,"%d-%dfinal.syn",myid, N);
			save_allocated_synapses(fname_syn);
			//save_sds(fname_syn);
			
		}
		//if (sec==1000) { // save synapses when we start to learn higher level area -- Pp23
		//	sprintf(fname_syn,"%d-%d_1000.syn",myid, N);
		//	save_allocated_synapses(fname_syn);
		//}
		
		// remove all spikes except those fired later than D ms ago
		k=N_firings; 
		while (k != B_firings) 
		{
			firings_time[k]=firings_time[k]-1000;
			if (--k<0)
				k = N_firings_max-1;
		}
		
#ifndef DEMO
		// write spikes
		//if(myid == 0 )remove("spikes.dat");
		char savefile[128]="spikes";
		char timestring[128];
		sprintf(timestring,"%d",simtime+1);
		strcat(savefile,timestring);
		strcat(savefile,".dat");
		if(myid == 0 )rename("spikes.datt", savefile);
#endif  
		
#ifdef Save_Dynamics_Every_Second

#if defined(PHASE_1_2_) || defined(PHASE_5_)
    /* Pattern + Rotation Training Settings */
    int save_interval = 1000;
#endif
#if defined(PHASE_3_) || defined(PHASE_4_) || defined(PHASE_6_)
    /* Blank-Out Test Settings */    
    int save_interval = 5000;
#endif
    if( (simtime+1)%save_interval == 0 )
      save_dynamics(fname_dyn);
#endif
		
		
	} // end of sec loop
	
	sprintf(fname_syn,"%d-%dfinal.syn",myid, N);
	save_allocated_synapses(fname_syn);

#ifdef LOG_GROUP
  //if(myid==0)
  fclose(filelog);
#endif

}



/*
 #include <signal.h>
 #include "singleton.h"
 #include "backtrace.h"
 #include "signal.h"
 
 namespace nsi {
 namespace spike {
 
 template <int Signal>
 class backtrace_signal
 : public signal<Signal>
 {
 typedef signal<Signal> signal_type;
 public:
 
 // Signal value
 enum value_type
 {
 signal_value = Signal
 };
 
 // Default constructor
 explicit backtrace_signal(bool handled = false)
 : signal_type(),
 handled_(handled)
 {
 
 }
 
 // Copy constructor
 backtrace_signal(const backtrace_signal& other)
 : signal_type(other),
 handled_(other.handled_)
 {
 
 }
 
 // Destructor
 virtual ~backtrace_signal()
 {
 
 }
 
 // Assignment operator
 backtrace_signal& operator=(const backtrace_signal& other)
 {
 handled_ = other.handled_;
 return signal_type::operator=(other);
 }
 
 // Create the handler
 bool create()
 {
 // This is required because we are not using fpermissive 
 // Clear action mask
 sigemptyset(&(signal_type::action_.sa_mask));
 signal_type::action_.sa_flags = SA_RESTART | SA_SIGINFO;
 signal_type::action_.sa_handler = 0;
 
 return signal_type::create();
 }
 
 // Operator which handles the signal
 virtual bool operator()(siginfo_t* info, void* ucontext)
 {
 ::nsi::spike::backtrace trace;
 trace.create();
 trace.print();
 
 return handled();
 }
 
 // Accessors
 bool handled() const
 {
 return handled_;
 }
 
 private:
 // If the signal is handled in the operator
 bool handled_;
 };
 
 void create_backtrace_handlers()
 {
 // Create our signal handlers
 backtrace_signal<SIGILL> ill;
 ill.create();
 backtrace_signal<SIGILL>::add(ill);
 
 backtrace_signal<SIGSEGV> segv;
 segv.create();
 backtrace_signal<SIGSEGV>::add(ill);
 
 backtrace_signal<SIGABRT> abrt;
 abrt.create();
 backtrace_signal<SIGABRT>::add(abrt);
 
 backtrace_signal<SIGTERM> term;
 term.create();
 backtrace_signal<SIGTERM>::add(term);
 
 backtrace_signal<SIGUSR2> usr2(true);
 usr2.create();
 backtrace_signal<SIGUSR2>::add(usr2);
 }
 
 } // namespace nsi
 } // namespace spike
 
 */

#ifndef STRINGIZE
#define STRINGIZE(x) (#x)
#endif // !STRINGIZE

#ifndef MAKE_STRING
#define MAKE_STRING(x) STRINGIZE(x)
#endif // !MAKE_STRING


#ifdef nomad
int main(int argc, char *argv[])
#else
void main(int argc, char *argv[])
#endif
{
	// Create our signal handlers
	//  ::nsi::spike::create_backtrace_handlers();
	
	chdir(MAKE_STRING(SIMDIR));
	try{
		// vvv MPI
		int numprocs;
		int  namelen;
		
		// Initialize MPI
		
		MPI_Init(&argc,&argv);
		MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
		if (numprocs<NP)
		{
			if( myid == watchid) fprintf (stderr, "The number of active processors, %d, must be >= NP=%d \n",numprocs,NP);fflush(stderr);
			exit (1);
		}
		
		// world rank
		int world_rank;
		MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
		
		// split
		int color = 'E'+'u'+'g'+'e'+'n'+'e';  // Other simulators must be different than this.
		cout << "************** first split color " << color << endl;
		cout << "Hello . . . " <<  myid << endl;
		// DEBUG: RGM
		int numprocs2, myrank;
		char processor_name[MPI_MAX_PROCESSOR_NAME];
		MPI_Comm_size(MPI_COMM_WORLD,&numprocs2);
		MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
		MPI_Get_processor_name(processor_name,&namelen);
		cout << "====STARK NNROBOT from " << processor_name << " - process " << myrank << " of " << numprocs2 << endl;
		// END DEBUG
		
		cout << "Before Split on process " << processor_name << " " << myrank << endl;
		MPI_Comm_split(MPI_COMM_WORLD, color, world_rank, &network_comm);
		cout << "After Split on process " << processor_name << " " << myrank << endl;		
		MPI_Comm_rank(network_comm, &myid);
		
		// split roots
#define MASTER 2
#define SLAVE 3
		
		color = myid == 0 ? MASTER : SLAVE;
		cout << "************** second split color " << color << endl;
		int key = 1;
		MPI_Comm_split(MPI_COMM_WORLD,
					   color,
					   key,  // I am the receiver of spikes; I will be procid 1 in this communicator; sender is 0;
					   &exchange_comm);
		
		MPI_Comm_rank(exchange_comm, &exchange_rank);
		
		MPI_Get_processor_name(processor_name,&namelen);
		
		if(myid == watchid) fprintf(stderr, "Starting\n");fflush(stderr);
		fprintf(stderr,"Process %d on %s\n", myid, processor_name);
		if(myid == watchid) fprintf(stderr,"processors: %d\n",numprocs);fflush(stderr);
		// ^^^ MPI
		
#ifdef DISPLAY_WINDOW
		int display_color;
		display_color = 101;
		int display_key = 1;
		MPI_Comm_split(MPI_COMM_WORLD, display_color, world_rank+1, &display_comm);
		MPI_Comm_rank(display_comm, &display_rank);
		int display_size;
		MPI_Comm_size(display_comm, &display_size);
		cerr << ">>>>>>>>>>>>>>> third split id " << myid << " : " << display_rank << " of " << display_size << endl;
		
#endif

// Check that only one(1) transmission option is defined !!!
#if defined(FIXED_TRANSMISSION) && ( defined(BITPACKED_TRANSMISSION) || defined(VARIABLE_TRANSMISSION) )
		cerr << "***** ERROR: Multiple MPI Transmission Options Defined! Exiting! *****" << endl;
		exit(0);
#endif
#if defined(BITPACKED_TRANSMISSION) && ( defined(FIXED_TRANSMISSION) || defined(VARIABLE_TRANSMISSION) )
		cerr << "***** ERROR: Multiple MPI Transmission Options Defined! Exiting! *****" << endl;
		exit(0);
#endif
#if defined(VARIABLE_TRANSMISSION) && ( defined(BITPACKED_TRANSMISSION) || defined(FIXED_TRANSMISSION) )
		cerr << "***** ERROR: Multiple MPI Transmission Options Defined! Exiting! *****" << endl;
		exit(0);
#endif
	
		// initialize GSL random generator
		gsl_rng_env_setup(); // GSL
		rand_generator = gsl_rng_alloc (randgen_type); // GSL
		
		
		//motor_world.init_robot(myid);
		srand(101);//RGM: ensure all procs have the same seed for DA value system coin flip
#ifdef EXP3
		::nsi::thinklink::scoped_ptr< ::nsi::thinklink::device > device;
		if (myid == 0)
		{  
			//must be cleaned up! RGM
            argc = 6; //from 6 to 9
            char msg[] = "t";
            argv[0] = msg;
            char msg1[] = "1";
            argv[1] = msg1;
            char msg2[] = "-a";
            argv[2] = msg2;
            char msg3[] = "10.11.2.253";//"apex-backpack01-wired"; //updated for wired connection IP
            argv[3] = msg3;
            char msg4[] = "-p";
            argv[4] = msg4;
            char msg5[] = "9675";
            argv[5] = msg5;
            //char msg6[] = "-t";
            //argv[6] = msg6;
            //char msg7[] = "3";
            //argv[7] = msg7;
            //char msg8[] = "-o";
            //argv[8] = msg8;
            
            ::examples::local_options opts(argc, argv);
            if (!opts.parse())
				printf("\n");
            
			std::cout << "*******************************************************************************" << std::endl;
            printf("Using socket device %s %u\n",opts.socket_address.c_str(), opts.socket_port);
            device.reset(new socket_device(opts.socket_address.c_str(), opts.socket_port));
            
			std::cout << "Opening Socket Device . . . " << std::endl;
            device->open();
			std::cout << "Socket Device is now open!" << std::endl;    
            g_motor_world::Instance();
            g_motor_world::Instance()->init_robot(myid, device.get());
			std::cout << "*******************************************************************************" << std::endl;
		}
		else
		{
			g_motor_world::Instance()->set_procid(myid);
		}
#endif
		
		
#ifdef EXP2
		motor_world.init_robot(myid);
#endif
		
		
		initialize();
		
#ifdef EXP3
		g_motor_world::Instance()->init(group_data);
#else
		motor_world.init(group_data);
#endif
		
		
		//run(start_sec+motor_world.get_sim_end_time_sec()); // run simulation for whatever seconds. 
		
#ifdef APEX3
    //initialize kbhit
    init_kbhit();
#endif

#ifdef EXP3
		run(g_motor_world::Instance()->get_sim_end_time_sec()); // run simulation for whatever seconds.
#else
		run(motor_world.get_sim_end_time_sec());
#endif

#ifdef APEX3
    //terminate kbhit
    kill_kbhit();
#endif
	
#ifdef EXP3
		g_motor_world::Instance()->reset();
#endif
		
		// destroy GSL random generator
		gsl_rng_free (rand_generator); // GSL
		
		MPI_Finalize();
		
		
	}
	catch (...){
		/*
		 nsi::spike::backtrace back;
		 back.create();
		 back.print(stderr);
		 */
	}
	
}
